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ABSTRACT 


We investigate an inverse imaging problem for TWI (Through-the-wall Imag¬ 
ing) using frequencies between 500 GHz and 1 THz. Starting from first principles, this 
thesis uses Maxwells equations to develop a model for the transmission Greens func¬ 
tion. This simplihed model is then used in a Lippman-Schwinger integral equation to 
predict the scattered held associated with interrogating THz waves. We investigate 
the effects of wave propagation through isotropic media, and present methods for cre¬ 
ating images from the scattered held. These methods are examined using simulated 
data. 
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I. INTRODUCTION 

Far better an approximate answer to the right question, which is often 
vague, than an exact answer to the wrong question, which can always be made 
precise. John W. Tukey, 1962 

X-ray computed tomography (CT) was created in the late 1960’s to early 
1970’s. Since that time, this imaging system has come to be recognized as one of the 
greatest achievements in radiology [Ref. 3]. From this point in time to the present, 
solving inverse problems and ill-posed problems has become an exciting and signihcant 
held for research. As in reverse engineering approaches, an “inverse problem” is 
simply one which attempts to attain information about a “source” from available 
information. 

Unlike a “direct problem” in which the descriptive parameters of the problem 
are known in advance of the solution, and the system of equations can be solved (sub¬ 
ject to initial and boundary conditions), an inverse problem poses a special challenge 
because the parameters are not known up-front. 

Suppose we want to measure the electromagnetic held at a given position 
caused by a signal x, and the held data d can be written as 

d = {Ax} (I.l) 

where A is an operator that describes the nature of the measurement system. 

Now let us entertain the idea that the object x G is mapped into the data 
space d G R™ via the operator A which would be a mapping of the object from the 
Hilbert space X into another Hilbert space Y. This problem is considered well-posed 
if it meets the conditions in Table I: 

Conditions 1 and 2 tell us that A is injective and spans all of Y. The require¬ 
ment of condition 3 is a necessary but not sufficient condition for the stability of the 
solution. The third condition is motivated by the fact that in all applications the data 
will be measured quantities. If the problem is well-posed, the error propagation is 
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1. A solution exists 

2. The solution is unique 

3. The inverse mapping is continuous 


Table I. Conditions for well-posed problems 

controlled by the condition number. For example, if there are n variables, there must 
be n independent equations [Ref. 30] in order to obtain a solution. Therefore the 
mapping of n variables by the system of equations to the data space can be written 
as Axn = d where A is a matrix dehning the system of equations and £c„ is a vector 
of n variables. By condition 1 and 2, A~^ exists. 

Now let’s say that there is a small error Ad in the measured data. Then the 
corresponding uncertaintyAcc^ for the object Ax^ obeys 


\AXr 


Xr 


< cond(A) 


l|Ad|| 

lldll 


The condition number is dehned as 


cond(A) = ||A|| ||A 


-ii 


where the norm 11A11 is 


I \Ax., 

= sup —j— 

Xn^O ||®n. 

If cond(A) is not large, then the problem is considered to be well-conditioned, 
and the solution will be stable with respect to small variations of data. If the condition 
number is large, we classify the problem as ill-conditioned. The difference between ill- 
posed and ill-conditioned is therefore somewhat vague. Hadamard dehnes a problem 
to be ill-posed if it does not satisfy all three conditions of Table I. An ill-posed 
problem is one where an inverse does not exist because d -|- Ad is outside the range 
of A. Generally, more than one image can map into the same data, or small changes 
in the data can result in large changes in the image. 
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Most inverse problems are ill-posed. Therefore, no unique answer exists. Let 
us consider the Fredholm integral equation, which we will discuss in depth later in this 
thesis. In real-world problems, the integral has a square-integrable (Hilbert-Schmidt) 
kernel and is of the form: 

k{x, s)f{s)ds = d{x) a < X < b 

If / has a small variation, say of the form A/ = esin(27ra;s), a; = 1,2,3,..., e G R 
and is small. Then the variation of data will be 

fb 

Ad{x) = e / k{x, s) sm{ 27 rcus)ds 

J a 

In Appendix A, we show that as a; ^ co then Ad ^ 0. Therefore, the ratio ^ 
heads to co as a; heads to oo. This problem is ill-posed, and the example demonstrates 
that Fredholm integral equations can be sensitive to high-frequency variations. 

Because of this behavior much research has been done in solving these types 
of problems, and many methods and techniques have been developed. This thesis 
will explore some of these techniques and create a simple electromagnetic scattering 
model for application to an inverse problem for high frequency imaging, appropriate 
to frequencies ranging from 10 GHz - 1 THz. 

A. BACKGROUND 

In recent years, terahertz technology has become an exciting field. Demands 
for further development have increased many fold for military, security, civilian and 
medical applications. As solid state physicists are coming up with equipment that can 
provide stable terahertz radiation, and as detectors are getting more sophisticated 
in measuring these data, scientists and engineers are now starting to focus on the 
terahertz region with excitement. This is largely due to the fact that: at these 
frequencies, there is great material penetration capability; unlike X-rays, the radiation 
is non-ionizing; unlike ultrasound, you can image without contact; and penetration 
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depth is generally greater than near-IR radiation [Ref. 8]. The most useful region of 
this electromagnetic spectrum is from 0.3 to 20 THz[Ref. 16]. 

B. SOME PROPERTIES IN THE GIGA-TERAHERTZ 
FREQUENCIES 

Frequencies of terahertz (a period of Ips) have some important properties. 
Most electrons in highly excited atomic Rydberg states orbit at this frequency. Col¬ 
lisions between gas-phase molecules at room temperature last about 1 ps. Super¬ 
conducting energy gaps are found at THz frequencies. Important collective modes 
of proteins are known to vibrate at THz frequencies [Ref. 16]. These are just a few 
examples of terahertz characteristics in the physical world. This region lies above 
the frequency range for most electronics but below the range of optical and infrared 
equipment; therefore, it has been difficult to explore until recent developments in 
solid state technology. 

However, there are limitations in propagating at terahertz frequencies because 
they are highly susceptible to atmospheric absorption. Only a small window within 
this region is practicable. 

The Terahertz Imaging Focal Plane Technology (TIFT) program at 
DARPA (Arlington, VA) dehnes the upper limit of this window to be 557 
GHz. 

There are also certain frequencies beyond this optimum band that can provide good 
propagation behavior. These windows are centered on 650GHz, 850GHz and at 
1.5THz. Anything above 1.5 THz results in absorption that is currently too high 
to gather any meaningful data[Ref. 8]. 

C. MOTIVATION 

As devices provide reliable radiation at THz frequencies, there is also a need 
for more sophisticated scattering models. There are many models already in place 
for other helds such as radar, and we have had much success in extracting superb 
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images in SAR, X-ray tomography, sonograms, etc. If these powerful methods can be 
carried over to the terahertz range, it may become both a rewarding and a worthwhile 
goal for development in the area. Recently, through-the-wall high frequency imaging 
has become a promising new technology that addresses many civilian and military 
problems. Much of the practical challenges he in the reliability of the system, its 
portability, and its ability to sustain an acceptable level of performance. 

One of the greatest challenges in imaging at these frequencies is that the system 
parameters are affected by the structured ambiguities in the walls, and this can cause 
difficulty in measurement. One has to account for diffraction and refraction effects, 
but the fact that we can image with relative success has shown that this area has a 
promising future. 

There are other factors that make terahertz radiation an attractive area for 
research. Waves at this frequency can penetrate non-metallic and non-polar media; 
moreover, explosives, chemical agents, and many different biological agents return a 
unique spectra. Furthermore, this radiation is safe for the operator of the device, and 
is harmless to the biological life that it may scan. 

Until the devices that generate terahertz signals can be made portable, prac¬ 
tical applications are still over the horizon. Academic institutions such as Rensse¬ 
laer Polytechnic Institute have demonstrated the reliabliity of these sources and are 
currently developing a more portable system for which identihcation, imaging, and 
various real world applications can be used. 

However, there are few other devices that are already in place that are small 
and reliable. Argonne National Laboratory (ANL) is currently investigating whether 
such penetration capability can be achieved at lower frequencies (in the 10 — 500 Ghz 
region). Such waves would be carried through a wave guide, and the return helds 
are measured via a waveguide. ANL demonstrated that these waves do have great 
penetrating capabilities, and they are able to form images. But imaging performance 
at these frequencies is an open question. 
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Though imaging can be done, most methods take advantage of imaging via 
bi-static transmission where the transmitter is on one end, and the receiver on its 
opposite end so that the image can be done via intensity strengths of held going 
through the media and the object of the image. However, imaging from the back- 
scattered held is more of a challenge because of the problems encountered not only in 
the perturbation caused by the object but also the interactions that occurs through 
the transmission medium. 

This thesis will emphasize the modeling requirements using Maxwell’s equa¬ 
tions and will develop inversion techniques for image reconstruction. We will hrst 
discuss the basic Inverse Problem. Then we will derive (from Maxwell’s equations) 
the theory for the simple model used in this thesis; derive the Green’s function for 
each region in space; establish the scattered held equation (perturbation of the held 
due to the object); discuss the various different approximations that can be done to 
establish the Fredholm equations for the scattered held; use the Born approxima¬ 
tion and parametize the equations; discuss the Singular Value Decomposition and 
the Least Square Method; establish the problems of error due to noise and discuss 
regularization; and hnally conduct some image reconstruction using simulated data. 
The hnal chapter will summarize the results of this thesis and provide some ideas for 
a more in-depth model for image reconstruction. 

D. PRELIMINARY DISCUSSION 
1. Material Effects 

For our model, we assume that the wall characteristics are known with a 
linear, homogeneous, and isotropic dielectric constant. There are certain ranges of 
Terahertz that are not very useful. For example, ranges above 5 to 15 THz frequencies 
have highly strong phonon absorption in virtually all crystalline materials [Ref. 11]. 
Innately metallic targets cause full terahertz reflection. Skin, due to the high water 
content, absorbs the terahertz rays, and most of the energy is dissipated in the hrst 100 
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microns [Ref. 12]. One must balance between the image object and the frequency that 
would provide the optimum penetration of the medium while providing the optimum 
reflection of the object. If we are looking for metallic material, it is an easier case 
because the focus would then be penetrating the medium. 

2. Atmospheric Effects 

Another important consideration is how feasible is it to transmit the electro¬ 
magnetic waves at these frequencies. There will be transmission losses in the held’s 
strength through the material and reflection of the target, but if atmospheric ab¬ 
sorption is too strong, one may just as well abandon the idea of imaging at these 
frequencies. Absorption and scattering in transmission occur largely due to the inter¬ 
actions of many different elements in the atmosphere. However, absorption is mainly 
caused by oxygen and water vapor: the greatest amount of absorption occurs when 
the frequencies match the natural resonances of oxygen and water vapor. As the 
frequencies deviate away from the natural resonant frequency of the elements, “win¬ 
dows” (frequencies at which absorption is minimal) are created. These are the areas 
in which we are especially interested. The Liebe model gives us a simple view of the 
attenuation and dispersive affects in the atmosphere [Ref. 14]. In Figure 1, we can see 
that frequencies between 0.7 THz to 1 THz have some great areas where attenuation 
would be small enough that we can neglect absorption effects. 

If we work at frequencies < 1 terahertz (such as 0.8 THz), we have almost 
no dispersive delay, and our attenuation is also quite small (< O.lOdB/km). Because 
atmospheric affects determine the type of devices and the methods of gathering data, 
we have to hrst determine which frequencies would give us: (1) the least attenuation; 
(2) the least dispersive effects; and (3) the most stable and common atmospheric 
conhguration. In Figure 2, we see that absorption is too high for any frequencies 
higher than 1 THz, and it will remain so until we get into the infrared and visible 
range — we will not be able to get any meaningful data in this region. Therefore, 
most of the research and work should be within this limit. 
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Figure 1. Attenuation a and delay P rates for three fog events {w = 
0.10,0.25, andQ.bg/vnP) added to a saturated air (RH=100%) sea level condition. This 
Figure is taken from [Ref. 14] 
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Figure 2. Atmospheric Attenuation from [Ref. 15] 

E. MODELING CRITERIA 

Now that the unexplored electromagnetic spectrum is starting to become acces¬ 
sible, the areas of research are growing without bounds. For behind-the-wall imaging, 
the interest for us is to develop a model that could predict held information using 
the Green’s function which, in turn, yields the system (PSF: Point Spread Function). 
The utility of having an analytical form for the PSF is that this is the kernel of the 
Fredholm Integral equation. Therefore, knowing the nature of the kernel which pro¬ 
vides the mapping behavior from one Hilbert Space to the other, we can then hnd 
the source (object) of the perturbed helds. Since we can measure the Helds, we can 
solve for the object x in equation I.l. 

To solve for the PSF, we will make the following assumptions: 

1. The waves are not dispersive: they each pass through the medium with speed 

dk ^ 

2. The medium is linear and isotropic. 

3. The medium is non-magnetic. 

4. The medium is non-conductive. 
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Note that the last two conditions allow us to simplify Maxwell’s equations. 
By making the transmission medium non-conductive, we allow the index of refraction 
to stay real and no damping will take place. This is a rather simple model, but one 
that can be adjusted in the future to implement this very real effect. Making the 
medium non-magnetic will make the slowing of the wave to be only dependent on the 
dielectric constant. And lastly to prevent the coupling of Maxwell’s equations, we are 
going to make the medium linear and isotropic. 
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II. 


THE INVERSE PROBLEM 


This chapter is intended to be a basic overview of the inverse problem that 
this thesis will cover. We will review the basic ideas of ill-posed problems and regular¬ 
ization methods for forming stable approximate solutions. We will conhne ourselves 
to linear equations with compact operators in Hilbert space. We will base our recon¬ 
struction algorithm design on the Singular Value Decomposition and the Tikhonov 
regularization method. There are many regularization methods for which the inverse 
problem for this imaging criteria could be done. The intention of this thesis would 
limit the scope of the most basic Tikhonov regularization method and the optimum 
methods for imaging would be achieved stochastically. For more in-depth discussion 
of ill-posed problems and different regularization techniques to enhance imaging, refer 
to Tikhonov and Arsenin [Ref. 26]. 


A. CONCEPTS OF ILL-POSEDNESS REVISITED 

In the last chapter, we saw that a well-posed problem must meet the three 
conditions listed in Table I. We now examine the data model equation in the form of 
a Fredholm integral equation: 

A{xi, s)f{s)ds = d{xi) a < s < b c < Xi < X 2 < ■ ■ ■ < Xn < d 



The interval [c,d] is partitioned into N parts for convenience to represent discrete 
measurements. This integral can be simplihed into the following matrix problem: 

f 1 

d{xi) = < lim A{xi,a + jAs)f{a + jAs)As > , where As = (b—a)/M (II.1) 

I A/f —»-nn • I 


j=i 


or 


d = lim Af 

M—>00 


where 
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d{xi) 


A(xi, a -|- As) 

y4(xi, a-|-2As) 

A{xi,a + MAs) 

d = 

d{X2) 

, A = 

A{x 2, a + As) 

y4(x2, a-|-2As) 

A{x2,a + MAs) 


d{XN) 


A{xn, a + As) 

A(a;Ar, a-|-2As) •• 

■ A{xn,cl + MAs) 


f{a + As) 
f{a + 2As) 

/(a + MAs) 


A will define the kernel which is to operate on / G X where X belongs to the 
Hilbert space. From the three conditions of Table I, we say that the mapping of A 
is therefore one-to-one and spans all of the vector space Y. Because we are dealing 
with compact operators [Ref. 27], Y will also be a member of the Hilbert Space. 
This tells us that the inverse A~^ must be continuous and bounded. A well behaved 
problem requires that / G X and the data d G Y are well controlled by the norms 
of the vectors in X and Y. But in most real life applications, we do not have this 
continuous inverse mapping, which suggests the contamination of the data space with 
elements not belonging to A.(X) C Y. This results in the instability of the solution 
for f which, as discussed in the previous chapter, is caused by the condition number. 
Colton [Ref. 27] proves that Af = d is improperly posed if / G X is not of hnite 
dimension. In everything we encounter in the real world, / is inhnite in dimension, 
and therefore, ill-posedness of the equation H.l is unavoidable. 

Straightforward solutions of ill-posed problems often result in non-physical 
answers. For this reason, one must be very careful in the discretization of equation 
H.l. One may think that making the operator A more rehned would make the 
solution more reliable, but the opposite is true [Ref. 27]. Therefore, techniques must 
be introduced to construct stable solutions. 


As 
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B. METHOD OF LEAST SQUARES 

A standard technique in mathematical and statistical modeling is to hnd a 
“least squares” curve ht to a set of data points. This technique was hrst developed 
independently by Carl Friedrich Gauss and A. M. Legendre. Whenever we measure 
data, we are going to get errors in measurements or inaccuracies. Therefore, we 
cannot expect the curve to perfectly ht through all the data points. Instead we want 
the best approximation in the sense that the sum of squares of the range between 
the data point and its corresponding curve is minimized. A least squares problem is 
generally formulated for overdetermined systems. For example, suppose x G M"' and 
b G M"* are related by the system of equations b = Hx. In general, we may or may 
not hnd a vector £c G M" to hnd b G M'”. We can, however, hnd a vector £c' whose 
transformation by H would result in the minimum distance between b and Hx'. 



H(X) 


Figure 3. 6 G while a:' G X and X = 
We form the residual: 


r = b — Hx' 


and we minimize the distance 


||r|| = \ \b — Hx'W 

From the geometry, we see that the minimum distance between b and Hx' has to be 
orthogonal to the plane. It can be shown by using Pythagorean Law that this will 
result in the least squares solution [Ref. 7]. From Figure 3, one can see that Hx' is 


13 






the projection of b onto the H[X). Due to orthogonality, r must he in the nullspace 
of H, where a null space is dehned as the set of all vectors a; G M"" such that for a 
given operator A, Ax = 0. In set-builder’s notation we have: 

Null{H) = {x : Hx = 0} (11.2) 

By the Fundamental Subspaces Theorem [Ref. 29] which says if A is an m x n 
matrix, then Null{A) = Range{A^)^ and Null{A^) = Range{A)^ and equation 
II.2, H^r = 0. Therefore, to solve the least squares problem Hx = b, we must solve 

r = b — Hx' 

H^r = H^{b - Hx') 

0 = H^b - H^Hx' 

H^b = H^Hx' 
x' = {H^H)-^H^b 

and is a projection matrix of b (otherwise known as the Moore Penrose 

Inverse). We will apply this method of model htting using regularization technique 
in the next section. 

C. BRIEF DISCUSSION OF REGULARIZATION METH¬ 
ODS 

In our previous discussion, we saw that the residual falls nicely into the 
nullspace of H^, but unfortunately in “real life” problems, this residual will not 
completely fall under the null space because it will be contaminated by some other 
element that we now include as noise. Methods for constructing a stable approximate 
solution of an ill-posed problem caused by the effects mentioned are called regular¬ 
ization methods. If error/noise is introduced into the equation d = Af such that 
d! = Af -|- n, then the perturbation would be: 

||<i' — d|| < ||n|| (II-3) 
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d will belong in the range of A, but d' may or maynot be in the range of A. From 
functional analysis, we introduce a linear compact operator R such that: 

lim RaAf = f a > 0 (II-4) 

This is called a regularization scheme for the operator A [Ref. 27]. a is considered 
the regularization parameter. For example, equation II.4 implies 

lim Rad A~^d 

As a goes to 0, Ra converges to A“b 

Colton [Ref. 27] points out that the regularization schemes for compact op¬ 
erators are not uniformly convergent. Because physical objects are generally infinite 
in dimension, RaA cannot converge to its norm as a —^ 0. For this reason, we must 
search for a good regularization parameter. The regularization scheme approximates 
/ of the data model Af = d with a regularized object f' satisfying 

f = Rad' (II.5) 


Including measurements, we have: 

f-f = Rad' - f 
f - f = Rad' - Rad + Rad - f 

f - f = Rad' - Rad RaAf - f 

f'-f = Ra id' -d) + RaAf -f (II.6) 

where d' — d = n, the (unknown) measurement noise. By the triangle inequality, we 
get 

II/' - /II < ||n|| \\Ra\\ + \\RaAf - /II (II.7) 

The error in object estimation depends on two terms of equation II.7 [Ref. 27]: The 
hrst term demonstrates how the error is propagated by the measurement limitation 
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and the noise of the system; the second term shows that error comes from the inverse 
mapping A~^ from specihc elements in Y mapped by A to X, and the linear operator 
R that maps all elements in Y to X. We desire that as the regularization parameter 
a goes to 0, so will the term \ \RaAf — /||. Unfortunately, this destabilizes the hrst 
term of eqn II.7. Since Ra —> A~^, ||J?q|| becomes larger (and thus the condition 
number). This behavior rapidly increases the reconstruction error. The hrst term 
in eqn II.7 therefore requires a larger regularization parameter for stabilization, and 
the second term requires the parameter to go to 0 for accuracy. Therefore, a method 
must be reached to optimize the right hand side of equation II.7. 

One way to accomplish this is to make a choice that the regularization param¬ 
eter a be dependent on the noise level, that is, a = a{n). From equation II.3, we see 
that d! ^ d as n ^ 0, and so 

Ra{n)d' A~^d = f 

There are many approaches to accomplish this behavior. In practice, a has 
to be determined via a priori or a posteriori criteria. In most of the cases for TWI 
(Through-the-Wall-Imaging), we are not certain what we are looking for, which makes 
a priori less practical than the a posteriori choice. For proof of concept, this the¬ 
sis will generate data with point scatterers (series of Dirac Delta functions), which 
we would then use the analytic form of the kernel to determine the point scatters 
themselves. We assume we have some prior knowledge of the model, which would 
not likely be the case in a “real world” type application. We will not cover a more 
detailed Bayesian statistical analysis of the model in this thesis, which is required for 
a posteriori criteria, but we leave it for future research to improve on the model. 

D. TIKHONOV REGULARIZATION 

We note from the previous section that as the regularization parameter goes 
to 0, we have an instability of ||i?a|| as a —> 0. This is primarily due to the fact 
that the noise contamination of the data is related to the (nonzero) eigenvalues of 
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Ra. As we try to get a more accurate solution, the first term of the error in equation 
II. 4 grows like (for proof refer to Colton [Ref. 27]). This is what some refer to 
as the penalty term where we punish the equality for a large norm. This form of 
regularization is due to Tikhonov [Ref. 26]. 

There are a few things we can do to minimize this effect. We seek to minimize 
equation II. 7 in the least square sense. 

= argmin{||n|p + \\RaAf - /|p} (II.8) 

all/ 

Using the fact that a is a function of noise, we can introduce this regularization 
parameter into equation II. 8. Consequently with this and equations II. 6 and II. 7, we 
can arrive at the Least Square solution for R^ via Newton’s Method [Ref. 27]. 

R^ = {al + A^A)-^A^ (II.9) 

Note that as a ^ 0, equation II. 9 gives us the Moore Penrose inverse, the 
simple least-squares solution. From SVD, we know that A can be written in terms of 
right singular vectors {V matrix) and left singular vectors {U matrix) with common 
singular values such that 

A = USV^ 

where /S' is a diagonal matrix with shared singular values (see Appendix C). Because 
U and V matrices are formed from eigenvectors of symmetric matrices of AA"^ and 
A'^A respectively, U and V are orthogonal matrices. 

A^ = {USV^f = VS^U^ 

Now we make the subsitution of A and A^ into equation II.9: 

= {aviv^ + vsu^us^v^y^vs^u^ 

Because U and V matrices are formed from eigenvectors of symmetric matrices of 
AA"^ and A'^A respectively, U and V are orthogonal matrices. Therefore, we get 
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the following: 

= {aviv^ + vs^v'^y^vs^u^ 

= [V{al + 5 ^)^^] VS^U^ 

= v{ai + s^yw^vs^u^ 



E. TRUNCATED SVD REGULARIZATION 

Another way to control the ill-posed behavior is to simply truncate the SVD 
representation: pick a small threshold value a and compare this value with a given 
singular value, which you will determine. Once this threshold is reached you replace 
the corresponding Da by 0. This method is called the truncation filter. There are 
many ways to determine where the threshold would he. In this thesis, we will look for 
where the steepest changes in singular value occurs and will truncated at that point, 
most of which would be trial and error. 
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Both methods—Truncation hltering and Tikhonov regularization—are effec¬ 


tive at reducing or minimizing the effects of noise. However, neither technique can 
independently offer a “best” threshold value a. 


19 



THIS PAGE INTENTIONALLY LEFT BLANK 


20 



III. WAVES IN MEDIA 

A. PRELIMINARY DISCUSSION 

In this chapter, we will set the foundations for applying Maxwell’s Equations 
to predict held perturbations in space located behind a transmission medium. We 
will work primarily with the time independent wave equations, particularly for the 
electric held. In free-space, the usual approach applies the free space Green’s function 
to the source term of the inhomogenous electromagnetic wave equation. This source 
is otherwise known as the scattering potential. 

If there is a source of perturbation that occurs in region III in Figure 4 and 
we have a receiver or measuring device in region I, generally the free-space Green’s 
function would work if region II is of the same medium of region I and III, that is 
free-space. However, if region II is a wall, the free space Green’s function does not 
satisfy the conditions determined by the boundaries of the wall. This will be ex¬ 
plained in the following sections where we break the spherical waves into half space 
expansions and apply the appropriate boundary conditions. We must also consider 
the effects of refraction and diffraction in the media. The Green’s function therefore 
needs to be solved either analytically or be measured. In real world situations, this 
function (also known as the System Transfer Function (STF)) cannot usually be de¬ 
termined analytically and is often measured. Such measurements are limited by the 
accuracy of the device. As discussed in chapter 2, noise contamination can result in 
catastrophic image reconstruction artifacts when null-space noise components con¬ 
taminate the data. Furthermore, the reconstruction becomes more unstable as the 
condition number increases. Robust and powerful techniques of regularization will be 
generally required. 

In the next chapter, we will create a simple model where we approximate this 
transfer function by imposing the appropriate boundary conditions. These results will 
then be applied (in the spatial frequency regime) to the development of an imaging 
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algorithm. 


B. PARAMETERS 

Electromagnetic waves through matter are governed by three properties of 
the material: the permittivity e, the permeability jj, and the conductivity a. In our 
model, we assume the material to be non-magnetic with constant permittivity and 
no conductivity, which means the index of refraction is going to be constant and 
real-valued. We will discuss some of the effects of the model and shortcomings due 

to this simplihcation. 

a. Maxwell’s Equations 

Maxwell’s equations for stationary media are of the following form[Ref. 

6 ]: 


V • E 


V B 


V X E 

V X B 


Ptotal Ptrue VP 

eo eo 

0 

dB 

f. dP ^ dE\ 

ho [jtTue + + V X M + ^0-^ j 


(111.1) 

(111.2) 

(111.3) 

(111.4) 


where ptotai is the total charge density composed of the charges that are free to move 
Ptrue and those charges that are bounded Pbounded = —V P, and P is the polarization 
of the medium, jtme is considered the true or the free current. And if we introduce 


I 


II 


III 


Figure 4. Region 
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the field vectors D and H by 


D = eoE + P 

H = —-M 
Mo 


then the Maxwell equations simplify to 


V • D Ptrue 


V B 


V X E 

V xH 


0 

dB 

dD 

3 true “I o. 


(IIL5) 

(III.6) 

(IIL7) 

(IIL8) 


1. Anisotropic Media 

For anisotropic media, the permittivity and/or permeability is not a scalar¬ 
valued parameter. Since most materials in the walls are non-magnetic, we can safely 
conclude that permeability will remain isotropic. However, most “real” problems 
may have anisotropic permittivity. Therefore D cannot be represented by a simple 
constant of proportionality to E. Instead a tensor equation is necessary so that [Ref. 
32] 

D = tensor{e) ■ E (III.9) 

where 

Bx ^xxEx T ^xyEy tixzEz 

E)y CyxE X “1“ ^yyBy “1“ ^yzE z 

Dz ^zxEx T (^zyEy -|- tzzEz 


and 


tensor{e) 


^xx ^xy ^xz 
^yx ^yy ^yz 
^zx ^zy ^zz 
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Since D may not be parallel to E, the permittivity tensor is a function of nine 
numbers. If we are considering a homogeneous media, then will be independent 
of position. Otherwise, they too can be varying. For most of the common cases, we 
expect the permittivity tensor to be symmetric [Ref. 32]. This means that only 6 
different constants are needed to dehne the 3x3 matrix. Furthermore, because the 
matrix is symmetric, its eigenvectors will be orthogonal. Therefore, if you choose the 
coordinate systems such that it is aligned to these eigenvectors, you can reduce the 
tensor to a diagonal tensor whose elements are its eigenvalues such that 


tensor'(e) 


400 

0 e'y 0 

0 0 e; 


(III.IO) 


This coordinate system is called the principal system. When the electric held lies 
along one of these principal axis (eigenvector), then D || £J,and thus the mathematical 
analysis can be simplihed. 

For simplicity, our model will neglect this effect (that exists in materials such 
as crystals). We will restrict ourselves to a linear media where the polarization of the 
media is directly proportional to the held. 


D = e^E + eoxE = eo(l + x)E = tQtrE = eE (III.11) 

where y is the electric susceptibility, cq is the permittivity of free space, is the 
relative permittivity, and e is the dielectric constant. 

2. Conductive Materials 

In a conductive media, which largely obeys Ohm’s Law, and assuming the 
material is linear, we write 

j = aE = -D (III.12) 

e 

By the equation of continuity 

V-j + ^ = 0 (III.13) 
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we see that 


= 0 

3 stationary 

SO that V • j stationary = 0 

Substituting equation 111.12 into this last result yields 

. = • 

3 ^ 0^ 3 stationary 

3 3 stationary 3 0^ " (III. 14) 

where Trei = “ is the characteristic or the relaxation time [Ref. 6]. 

Taking the curl of equation III.7, we obtain the following 

d 

V x{V xE) = -^(V X B) (III.15) 

We also restrict our model to be non-magnetic. Most wall materials are indeed non¬ 
magnetic, and we are safe to employ this approximation. Therefore, the relative 
permeability ~ 1 and so the material permeability is approximately the magnetic 
constant in vacuum jj, « yUo- This allows us to write B = HoH. With Equation III.8, 
Equation III. 15 then can be simplihed to 


V- j 


3 + 


dD\ 

) 

dD 


dt 


V(V • E) - V^E = -yuo^ [jtrue + j (111.16) 

For most cases, we will be dealing with charge free regions. Hence \/ ■ E = 0. 
Now for conductive material, we may substitute Equation III.12 into III.16 for further 
simplihcation. We arrive at the following with the given assumptions: 

d'^E dE 

Harmonic solutions are of the form: 
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E{x, y, z, t) = E{x, y, z)^e 


From this point, we shall work with time independent wave equations, and 
will refer all vectors in the frequency domain, and Equation 111.17 becomes: 


“ 1 “ iJiQtuj^E -|- iyodOjE — 0 

V^E + = 0 


(111.18) 

(111.19) 


Since e = eoe^, c = 



and n 


we can simplify the above equation to 


y^E+ ( 

^ icr\ 

0 

1 H-I 

^oj^E 

\ 

V 


V^E + 

(i + A 

V eo; , 

) n^k^E 


(III.20) 


where uj/k = c. With the relaxation time derived earlier T^ei = e/a, we get 


V^E+ 1 + 


'^rel^ 


n^k'^E = 0 


(III.21) 


From Equation III. 18 we define a complex wave number k = k + in so that 


P = yoeuj^ + iyoduj 


(III.22) 


and 


k 

K 




1/2 


+ 1 



Note that for all pure metals, the relaxation time is « 10“^"^ [Ref. 6], which is on the 
same order as the collision time of the electrons. At the frequency that we propose 
to use (500GHz to ITHz), the relaxation time will be of the order « 10^ if the 
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material is fully conductive. This means that the diffusion part of the equation would 
dominate and therefore the attenuation would be quite high. The penetration depth 
governed hy d=l/n will be quite short, and the waves will not be able to penetrate. 
Usually for a good dielectric the relaxation time is of the order of days[Ref. 33], 
and hence the imaginary term of Equation III.21 vanishes, and thus there will be no 
damping in the media. For such frequencies, we maybe able to generalize to certain 
conductive materials with larger relaxation times, but such an analysis would require 
a large database of material effects and is beyond the scope of the thesis. For this 
reason, the model also assumes that the walls are non-conductive and we will deal 
with just the wave equation portion of Equation III.17 

3. Frequency Dependence and the Dispersion Relation 

Dispersive media have values of yu, e, or a that are functions of frequency. It 
can be shown that for a dielectric e as a function of frequency is well-mo deled by [Ref. 
33]: 


e{uj) 

e'{uj) 


— ie\io) 

{io^ - 
eouj^uja 


where a = 1/t and r is the collision times of the electrons, uq is the resonant fre¬ 


quency. 


Lul = - plasma frequency, 

^ Com 

N are the number of dipoles per unit volume, and m is the mass of the electron. The 
real part d{uj) dehnes the refractive index n(cu) = ^/e' {uj)/eo, and the imaginary part 
e"{uj) is the loss tangent of the material tan0(a;) = e"/e'{uj), and it is related to the 
attenuation constant (absorption coefficient) of the electromagnetic wave [Ref. 33]. 
In Figure 5, we see that around the resonant frequency cuq, the dielectric behaves in 
a strange manner and the material absorption becomes quite high. Therefore, one 
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has to be careful that the radiation frequency we transmit would not be in resonance 
with the dielectric, and further material analysis might be required. 



Figure 5. Real and Imaginary parts of dielectric constant. Figure taken from [Ref. 
33] 


For metals, charges are unbounded, so ujq = 0, and from Ohm’s Law, we hnd 
[Ref. 33] 

a{uj) = (III.23) 

In “real” phenomena, we know that from optics that n is a function of frequency 
(See Figure 6). Therefore, it can be seen that different wavelengths of light refract 
differently. This is largely due to dispersion dw/dk ^ constant. Because waves of 
different frequency travel at different speeds in a dispersive medium, one would have 
to find methods and ways to isolate these frequencies for application. The waves as 
a whole or the “envelope” of the waves travel at group velocity: 

dio 

Because the velocity depends on the frequency and uj = u!{k), the field is of 
the form: 


J —oo J —oo 

The amplitude A[k) is related to the initial condition of '0(x, 0) by: 

poo 

^(x,0)= / A(/c)e*^"dk 


(III.24) 


(III.25) 
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Analysis of the system typically requires a numerical approach (such as the the 
Method of Stationary Phase [Ref. 5]). The motivation for this is that we know that 
for all wave numbers k, wave component disperse for any initial conditions. All the 
waves will travel with the group velocity, and all of the energy maybe concentrated 
in one wave number [Ref. 5]. For this reason, using the Method of Stationary phase 
can provide us a very good prediction of the system. These are all physical realities 
we will have to face in dispersive mediums, for which case a large database of various 
different material will be required. It is beyond the scope of this thesis to analyze the 
effects of the waves for each different media. For this reason, this thesis will set the 
criteria that we are working with non-dispersive mediums. 

C. FOUNDATIONS OF THE SCALAR THEORY 

From our Maxwell equations listed in Equations (III.5), (III.6), (III.7), (III.8), 
we get our wave equations by taking the curl of Equation III.7 and Equation III.8. 
Since the vector wave equation is obeyed by both E and H, an identical scalar wave 
equation is obeyed by all components of those vectors provided we are in a cartesian 




Figure 6. Index of Refraction vs Wavelength. This Figure is taken from [Ref. 31] 
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coordinate system. Thus, the component will, for example, obey the equation: 


= 0 

Therefore, it is possible to summarize the behavior of all components of E and H 
through scalar equations of the form: 

+ = 0 (111.26) 

If we take from Equations (HI.7) and (III.8) and compare component by com¬ 
ponent, we will get 6 equations, and if we orient the coordinate system so that one 
component vanishes, then the equations will simplify to two equations with the E 
and B held. These equations can be uncoupled and reduced to a single scalar wave 
equation, which describes either a Transverse Magnetic Polarization (TM) or Trans¬ 
verse Electric Polarization (TE). When, all the medium parameters e, a, and jj, are 
independent of one of the space dimensions, then the corresponding electromagnetic 
held components are also independent of that same dimension. Then one can write 
the vector Maxwell equations as a superposition of the TE and TM helds. 

If there is variation of permitivity as in the case of anisotropic media, the wave 
equation becomes 


V{E-V{\ne))+V^E + n\^E = 0 (III.27) 

In this case, the extra term couples the x, y and z components. The scalar represen¬ 
tation would fail in the conductor case where the boundary conditions are that there 
is a jump of a in the electric held perpendicular to the surface, and a jump with the 
magnitude of the free current in the magnetic held parallel to the surface. In these 
cases, the components of the Electric and Magnetic helds will be further restricted, 
and would require a more robust model such as the Stratton-Chu model for vector 
wave equations [Ref. 27]. 
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If the waves are spherical in nature, one can still formulate the problem as 
a scalar wave in terms of the Debeye potentials [Ref. 6]. In either case, our model 
criteria gives us the flexibility for using the scalar model. Although there are faults in 
using the scalar model, our model criteria allow us to decouple the Maxwell equations 
and reasonably capture the essential behavior of the wave propagation. 

Furthermore, we have to impose the conditions at inhnity. In the vector case, 
the electric and magnetic helds must go to zero at inhnity, and thus must satisfy the 
Silver-Miiller conditons: 


lim [H XX — rE) = 0 (III.28) 

r—^oo 

lim {E X X rH) = 0 (III.29) 

r—»-oo 

where r = \x\ and where the limit is assumed to hold uniformly in all directions. 
The corresponding radiation condition for Cartesian components is the Sommerfeld 
radiation condition 

lim r(—- ikE) = 0 (III.30) 

1 -»^CXD QT 

For more complete discussion and proof, refer to [Ref. 27]. 

For the reasons mentioned above, we shall attack the problem for the scalar 
case with the Sommerfeld conditions imposed at inhnity. 
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IV. THE INHOMOGENEOUS WAVE 
EQUATION/TRANSEER EUNCTION 


To develop the transfer function, consider a field that emanates towards an 
infinite flat slab (Figure 7). We take this wave to be directed in the positive 2 ; 
direction, and y axis is positive downwards, and x is positive outwards. The isotropic 
media has an index of refraction that is independent of frequency, and the location 
of the source is at Vq = —zqU. 

For a held with source s(r,t), the wave equation becomes 

where £ is in the time domain for E. The source function can be written in terms of 
the Fourier integral: 



Figure 7. Point Source emanating waves through a constant dielectric medium (In- 
hnitely large slab) 
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sir,t)= / 5'(/,r)e*2-/Mf 


Thus the time independent wave equations reduces to 


+ n^eE = S{f,r) (IV.2) 

We can construct the solution of Equation IV.2 by superposition of unit point solu¬ 
tions corresponding to a source at the point tq = —zok, and we placed this source on 
the z axis for reasons to be discussed below. 

The Green’s function due to a source at tq will satisfy equation: 


(V^ k‘^n^{z))G{x, Vo) = S(x + Vq) 

where fc = 27r/A, A is the wavelength in free-space 

S{x Vo) = S{x)S{y)S{z -I- zo) 


and 


1 for 
n{z) = •^ n for 
1 for 

Therefore, the particular solution for E is 


z <0 
0 < z < L 
L < z 


(IV.3) 


(IV.4) 


(IV.5) 


E=[ G(a;,ro)S(/,r)d’r(i 
4 Vo 

and, in the time domain, we have 


£{r,t) 




(IV.6) 


(IV.7) 


Since the media is non-conductive, the wave number remains real valued. 
Therefore, there will be no dampening term, and the media will become a wave guide 
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composed of TE and TM waves. Using Bracewell’s notation, we take the Fourier 
transform of Equation IV. 3 in the x and y coordinate to obtain 

/ OO noo 

/ G{x, y, z, 

•OO J —OO 

/ OO /*oo 

/ G(/., /„ 

•OO J —OO 

This transform will reduce the three dimensional Laplacian to a single ODE. We 
are allowed to make this substitution because there is no slowing or changing of the 
media in the x and y direction: the media remains constant in these directions, and 
the index of refraction only changes in the z direction. Application of Equation IV.3 
yields 


^2 

[-{2'xf^f-{2'xfyf + —]G{f^Jy,z,rQ) + k‘^n^{z)G{f^Jy,z,ro) = S(z+Zo) (IV.IO) 


Dehne 


F = 27r fj + 27r fyj 

so that |Fp = {271 fxf + {271 fyY. 

The wave number is related as: 



= (27r/^)^ + {277 fyf + {277 


(IV.ll) 


(IV.12) 


Consequently, the Bracewell notation chooses the spatial frequency to obey = 
l/Ax ,fy = l/Xy,fz = l/Xz and we write G = G{f^, fy, z, Vq) 

When z ^ —zq, Equation IV.IO is the homogeneous wave equation with solu¬ 
tion 


G 


i^k^n^{z)-\F\^z ^ ^ 

N. y ' S. 

^ v >— 

right going wave 


■ix/ k‘^n‘^ {z) — \F \ ‘^ z 

V ^ 

left going wave 


(IV.13) 
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This solution must satisfy the boundary conditions, which will determine the coef- 
hcients of each component wave as they propagate through the medium. The wave 
and its derivatives must be continuous at any boundaries. Consequently, the following 
conditions must be true (see Figure 8): 


Y 



K 

6™ 

a 

(0,0,-z_0) 

d" 


A “ 


B 

G 

G 

G 


L 



Figure 8. Cross sectional view of the wavefront 


- 1 /^m 1 

12=0 — ^ |2=0 

OG- 

dG^\ 

dz 

II 

0 

Q:) i 

\z=L = G^\z=L 

dG^ 

dz 

dG^ 

z=L dz 


(IV. 14) 


(IV.15) 




— r^- 

Z=-Z(, 


dG- 


lz=-Z0 


dz 


dG^ 


Z— — Z0 


dz 


= 1 


(IV. 16) 


Z^-ZQ 


Imposing the Sommerfeld conditions and using the fact that the wave is outgo¬ 


ing from the source, the coefficient of the right going component (for 2 ; < —zq), and 
the coefficient of the left going component for [z > L) must vanish (see Figure 8). 
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The waves are traveling to the left when z < —zq and are traveling to the right when 
z > L. The boundary conditions are such that the Green’s function is continuous 
with continuous derivatives — except at the source. Since the waves are continuous, 
the Green’s function has to be continuous. But its hrst derivative is not very obvious. 
Let us integrate Equation IV. 10 both sides. 

G{fx,fy,z,ro)dz+ J (k^n^(z))-(27rfx)^-(27rfy)^)G(fx,fy,z,ro)dz = J h(z+zo)dz 

(IV.17) 

The right hand side will become the Heaviside Function. The left hand side would be 
composed of the hrst derivative of the Green’s function and an integral of the Green’s 
function. 

r 

fy, z, ro) + J {k‘^n^{z) - (27r/^)^ - (27r^)2)G'(/„, fy, z, ro)dz = H(z + zq) 

(IV.18) 

Since, the Green’s function is continuous, the middle term will vanish at any 
given boundary. However, the hrst derivatives will become continuous as long the 
boundaries are z < —Zo or z > —zq. As the hrst derivative crosses the source (where 
the Heaviside Function transitions from 0 —> 1, the difference across will have a jump 
of 1). And imposing these boundary conditions specihed above, we arrive to the 
following: 

Let 



At z = 0 


tit = n{z < 0) 

Tijn = n(0 < z < L) 
ub = n{L < z) 


a- + fc- = a™ + 


(IV.19) 
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sjk^nl - \F\^a- - - \F\'^ = 

|F|2 - b^^k‘^nl^-\F\^ 


1 1 


a 


1 1 



■\/ Nx —\/ Nb 


b- 


V^m —\/ Nm 




Ml M2 


where Nt = k‘^n^ — \F\'^ and Nm = k'^n^ — \F\'^. 

At z = L, we require 

^mgi\/fc2n^-|F|2L _|_ |^7n^-i^yk‘^n‘^-\F\‘^L _ k'^n%-\F\'^L 


(IV.20) 


(IV.21) 


(IV.22) 




(IV.23) 


NfnL g i'J NmL 


a™ 





M3 


(IV.24) 


^iy/NsL g—iV NbL 


’ ’ 

g-ivW^L 


0 


M4 


where Nm = k^rim — \F\'^ and Nb = k^n% — l-Fp 
And at z = —zq, we require 


y>^iy/k'^n^-\F\‘^zo _ Q-g-i\/ k^n'^-\F\'^zo _|_ ^i\/k^n'^-\F\‘^ zq 


(IV.25) 


Now we apply the jump condition (see Equation (IV. 16)) to obtain 
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(IV.26) 


a 





^-iy/NTZQ ^i^Ntzo 


a 



b- 


M5 


(IV.27) 


^i^Ntzo ^-iyJNxzo 

'-V- 

M6 

where A^t = — \F\'^. We note that all the matrices are invertible, and therefore, 

we can have a closed form solution for the coefficients of the Green’s function. 

To solve for these coefficients, we hrst simplify our results to read 



Ml 

a 

= M2 

o'" 


b- 



M3 

a™ 

= M4 

' ' 




0 



a 


1 

1_ 


Ms 


= M, 


+ 


b- 


0 



(IV.28) 


(IV.29) 

(IV.30) 


where the reflection and transmission amplitudes can be expressed in terms of a , 
the amplitude of the incident wave: 




B 

a 

= Mf ^MaMg-^Mi 

a 

b- 

'-V-' 

B 

0 


(IV.31) 
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or 


B 


a = -Bi.io 
b- = -®“ - 


B 


b = -02,10 

D d 

= 


B 


1,1 


B 


1,1 


After much algebra, we obtain the reflectance: 


(IV.32) 

(IV.33) 


^_b- _ 

a- (y]v^+^)(y]Vj:+^) + (-^+y]v;;)(y]V^-^)e* 2 ./]v;;:L 

(IV.34) 

The Reflection coefficient then is .R = r * r. 

Similarly, the transmittance is 


^ ~ (^/At + ^MnKVN^ + ^/]v;;)e-v^^ + 

(IV.35) 


The Transmission is T = t *t. 

We can also solve for the other coefficients. Inside the medium, we have 


a™ 

= M3-^M4 

' ' 


c 

0 


(IV.36) 


a™ = Ci.ia^ = 


C*!,! - 

Bi.i 


fo- = C2,ia^ = 


C2,l _ 

-°i,i 


(IV.37) 


a 


m 


a 


_ + _ 

+ y/KuKy/T^ + y/Kn) + {-y/N^ + 

(IV.38) 


a~ 


__ 

(vAVb + y/KnXy/Nr + y/Kn) + {-y/Nr + VKn){^ - 

(IV.39) 
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Finally, at the source, we obtain 


0 






0 


Writing X = Mg ^M 2 M^ ^M 4 and 



(IV.40) 


yields 


S 


0 


1 

e*\/W^o(i+/yy) 

1 


1 

e-*\/M^o(i+/;vy) 


(IV.41) 


b' = = (W,i/Sia)a- - (IV.42) 

and 

0 = X2,ia^ - 52,1 (IV.43) 

We can simplify these results somewhat. Outside the medium, we assume 
= V^B = \fN ■ We are interested in the amplitudes of the transmitted wave 
(a^), the incident wave (a~) and the reflected wave (b~). 

In the next chapter, we will show that we can consider the Green’s function 
to be shift-invariant so that G(x, Vq) = G(x — Xo,y — yo, z, zq) = G(i — ro). 

^i^/NzQ 



{N - N^) + {N^ - 


(IV.45) 


= t * a = 




i\fNzQ 


(IV.46) 


(v/iv + - {^/N^ - 2y/N 

We note that for 2 ; < —zq, is free space. Since the free space Green’s function is 
a solution to the Helmholtz equation and satishes the radiation condition at inhnity. 
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it must be the solution in this region. We will show this for —zq < z < 0 in the next 
chapter. For the regions of interest, we have the following: 

Zo < z < 0 right-going 
left-going 

(IV.47) 


Gifx,fy,z,ro) = 


2^/]V 

(\/v;7-i-A/iv)2-(^/]v;;;-^)2e»2yj^r 2^/]v ^ 


42-|F|2 



GifxJy,z,ro) = right going wave (IV.51) 

G{fxJy,z,ro) = left going wave (IV.52) 
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V. 


THE ANGULAR SPECTRUM 


A. WEYL’S INTEGRAL AND SPHERICAL WAVES 

In this chapter, we will show that the right-going wave of Equation IV.47 is 
the planar expansion of the spherical wave by introducing the Weyl’s integral. We 
will also discuss the properties of the Inverse Transform of the Green’s function in the 
spatial domain, and demonstrate the shift-invariant form of the Green’s Function. 

What makes the problem very difficult is the reflection and refraction of a 
spherical wave on a planar media because the wave has spherical symmetry and the 
boundary is a plane. For this reason, we want a way to expand the spherical waves 
in terms of planes. In equation IV.47, we will hrst take the case for the right-going 
wave. The Inverse Fourier Transform dehned as: 


G{x, y, z, ro) = / / G{f^, fy, z, ro)e*""(^^"+^-^)d4dfy 

J —oo J —oo 

H. Weyl’s expansion of the spherical waves into planar waves is [Ref. 24]: 


jk\r-r I poo POO ^ 


\r — r'\ 2 tt 


' —OO J —OO z 




where r = xi + yj + zk, r' = x'l + y'j + z'k and 


(V.l) 


(V.2) 


s, = +Jl- si- si when + sj < 1 


Sz = +iy sl + si — 1 when si + s^ > 1 
Therefore, the Inverse Fourier Transform for the right going wave becomes 


G{x,y,z,rQ) = 


POO ^i^yk'^-{27^fx)^-{27Tfy)^{z-\-zo) 

1-00 2ky^l - (27rA/A:)2 - {2'Kfy/ky 


(V.3) 


Under the change of variables: 
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_ 


dSr = 


„ _ 27 r/y 

k 


Zvrdfx 1 _ 27rdfy 

k “ k 


Equation V.3 becomes: 


noo i>oo k^-k"^sl-k"^ sl{z+zo) ^i[ks:cX+ksyy) f ]^ \ f 

G{sy„Sy,z,ro)= I I - 


' —OO j —OO 


2fcy/l -si -si 

^ /-OO /-OO ^i{k^l-s%-s^{z+zo)+ks3;X+ksyy) 


Svr^ 


' —OO J —OO 


\/l - 4 - s 


dSxdSy 


and, substituting equation V.4 into equation V.2, we arrive at the following: 


G{x,y,z,rQ) = - 


I gifc|r-ro| 


iAir \r — ro\ 

Operating at far field, this Green’s function would then be simplihed to 


jk\r-ro\ 


^ikr 


^-ikf-VQ 


i47r|r —ro| i47rr 

where f is the unit vector pointing in the r direction [Ref. 25]. 


(V.4) 


(V.5) 


(V.6) 


B. GREEN’S FUNCTION IN MEDIA 

We make a coordinate transformation of Equation IV.47 and Equation IV.48 
for reasons to be explained in the following sections. Use the relationship in equation 
IV. 12, and substitute the spatial frequencies in for the wave number such that we get 
the following (see Figure 9): 


27r/a; = fc sin 0 cos (f) ‘^T^fy = k sin 9 sin cj) (V.7) 

271= k cos 9 

In equation IV.47, we notice that ^* 2 ^° is the planar expansion 

of the spherical wave. We use this observation to see some physical similarities of 
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Figure 9. Coordinate Axis 


the other waves. The reflection of the wave in Eqnation IV.47 can be written as 
And so the returning wave is also spherical in nature: 


R = 


{N - NJ + (Nm - 

(^ a ;; + VNy - ( va ;; - 


(V.8) 


a/aJ^ = - (27rA)2 - {2'KfyY '/N = _ (27rA)2 - (27r/y)2 (V.9) 

From equation V.7, we rewrite eqnation V.9 as 

\/Nm = k\J — sin^ 9 \/A = k\/l — sin^ 9 = kcosO (A.IO) 

Substituting equation V.IO into V.8, we get the following: 

(1 - n^) + (n^ - ly^kLy/n^-sin^e 

ihyuJ — , 

{yj'n? — sin^ 0 + a/ 1 — sin^ 0)^ — (-\/— sin^ 9 — a/I — sin^ gy^i 2 kLy/n'^-sin^ e 

(V.ll) 

(1 - n^) + (n^ - i)ei2fcLVn2-siiP0 

( a/tt,^ — sin^ 9 + cos 6)^ — [\/ rfi — sin^ Q — cos Q^e^'^kL^/n'^-sm'^e 
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The Inverse Fourier Transform defined in equation V.l, with the following 


change of variables: 


fx 


— sin 9 cos (j) 

ZTT 


fy 


— sin 6 sin 6 
27r ^ 


so that 


dU 

dJjL 


ae 

a4> 

_ 

dfy 



ae 

a4> 



^ cos 9 cos (f) — ^ sin 0 sin < 
^ cos 6* sin ^ sin 9 cos 4- 


The spatial reflected Green’s function becomes 


— 1 cos 6* sin 6* 

ZTT 


G{x,y,z,ro)= / / G(6/, ] cos9sm9d9d(j) (V.12) 

' ' ' ZTT ' 


where G(6',z,ro) = -^(^) ""tklZe ^ 

It is easy to see that, under the coordinate transformation, Equation IV.47 
can be made into Equation V.12, but the limits of integration are not very obvious. 
Because we are summing the contributions of all the waves from various propagation 
directions, we can see that (j) varies from 0 to 27r. The projection of the propagation 
vector onto the x — y plane is symmetric for any given angle 9 normal to the me¬ 
dia. However, the 9 coordinate requires more careful analysis. In its cartesian form 
(Equation V.l), the Inverse Fourier transform suggests that there will be an area of 
integration where and fy, the projected components of the propagation vector onto 
the X — y plane, will be large enough that the argument of the exponential of z in the 
frequency domain of the Green’s function is no longer complex-valued but instead, 
becomes exponentially decaying in the z direction. In the following section, we shall 
discuss this effect, and discuss how 9 will have to be adjusted to include this effect. 


C. ANGLE OF REFLECTION/INCIDENCE 

Because the Green’s functions are continuous at the boundaries, the total field 
at the interface between any two media must obey incident wave + reflected wave = 
transmitted wave for all the regions in space where there is a change in the index of 
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refraction. Since this requirement must hold true for all time, it must be frequency 
independent. Therefore 


oj{z < 0) = ci;(0 < z < L) = cu[L < z) (V.13) 

Now, for the reflection at a surface, this requirement becomes 


c(^Z < Q'jkright—gcnng c(z < 0)fc;e/t— 


going 


k 


right—going 


= k 


left—going 


(V.14) 


Because the waves must be continuous at the interface, the projection of the prop¬ 
agation vector onto the x — y plane will be the same at the boundaries(z = 0 and 
z = L). Therefore 


kx = ‘2'Kfx, ky = 2'xfy, and fcsin^ 

must the same at z = 0 and z = L for all the waves (incident, reflected, transmitted). 
From Equation IV.47, a single wave propagating to the right from the source can be 
written as 

G'.Weni = G(/., fy, Z, 


„i\fN 20 , -j=j— 

2^/N 

and the reflected wave (the left going wave) is written as 


(V.16) 


G 


reflected 


{N - Nm) + - N)e 


{VKn + y/Ny - - y/Nyei^VN^L 2y/N 


i2\/NmL iVNzg ^-=— 

_ Q-W k'^-\F\^z^i{2'Kf^x+2-Kfyy) 


(V.17) 


Comparing the exponents, we see that only the 2 ;-component of the propagation vector 


is different by a negative sign with the same magnitude. And seeing that the k must 
be equal, we conclude 

VP - |F|2 = yp _ |^2| 

'---' vl--/ 

incident reflected 
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k cos 9inc = k cos 9ref (V. 19) 

9inc = 9 ref (V.20) 

This is the familiar Law of Reflection. 

D. SNELL’S LAW 

For transmission, we can verify Snell’s Law in the same manner: The frequen¬ 
cies and the projection of the propagation vector at the x — y plane need to be the 
same due to the boundary conditions discussed earlier. This requirement means: 


uj{z < 0) = C(;(0 < z < L) 


ckirtr — k 


n 


trans 


nkjrtr = k 


trans 


If we measure Otrans from the normal (see Figure 10), we can see from the 
geometry that 


Sm Ofj-ans ^mc Sm 0inc jTlky 

Ti sm Oifans sm OiYic 


(V.21) 


This is Snell’s Law. Equation V.21 suggests that there is a limit to the angle for 
transmission because the sine function has a range between —1 and 1. But if we need 
to integrate further into inhnity, we cannot be restricted to real angles. Stratton [Ref. 
34] introduces the complex angle of incidence through which the sine functions will 
become hyperbolic functions. This approach can be used to account for the evanescing 
waves, and is implicit in the Inverse Fourier Transform of Equations (IV.47), (IV.48), 
(IV.51),(IV.52). In the Cartesian form, we can see that fz ioo as A ^ Too or 
fy Too. 

To see this in the spherical transformation, we examine at the 2 ;-component of 
Equation V.7: cos 6* = (27r/^) ^ combination of /J and fy is 

larger than k^, then 
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Figure 10. Snell’s Law 


cos 6* 


iy/{27rf^y + (27r/^)^ 
k 




ICO 


(V.22) 


as either or fy approaches to inhnity. From the geometry in Figure 9, 6 must 
go from 0 to n/2, but which branch of 6 should be used in the complex plane? A 
complex angle of incidence indicates that 9 should be of the form 6 = a + iP, and so 
by trig identity: 


cos(6*) = cos(q: + iP) = cos(q:) cos(ip) — sin(Q:) sin(i/?) (V.23) 

= cos(q:) cosh(/?) — ■jsinh(/?) sin(Q:) (V.24) 

Our allowable angles 0 < 7r/2 and, in order for Equation V.22 to be true, P 
must go towards —oo since sinh(/3) is an odd function. Therefore, we choose a path 
for which 9 = 0 ^ f — too. 


E. EVANESCING WAVES 

We consider the region —zq < z < 0. From Equations (IV.47) and (IV.51), we 
consider the reflection and the transmission for a single wave front. 


G(x,y,z,ro) = A(A, 

G(x,y,z,ro) = A'(A, 

where A{fx,fy) and A'{fx,fy) can be complex or real-valued amplitudes. 


(V.25) 

(V.26) 
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Consider the reflected wave: the propagation vector of the reflected wave must 
be parallel to the x—y plane before the wave can start to evanesce. Once this condition 
exists, we can see that the exponent in the z direction becomes real, and the only 
imaginary portion is the argument 27rfxX + ^irfyy = F ■ xi + yj. Therefore, the wave 
propagation is parallel to the media before it can start to evanesce. But from a larger 
index of refraction to air (which is the case at 0 < z < L) where the rays bend away 
from the normal, one can see from Equations (IV.48) and (IV.51) that there is an 
angle of incidence will result in a transmission angle that is parallel to the surface. 
This is the condition for total internal reflection, and the angle is the critical angle. 
From this view point as the 2 ; component of the incident wave inside the medium 
decreases, the angle of incidence increases, and thus the transmitted wave will start 
to evanesce. As either or fy increases further, the conditions are met where the 
reflected waves will also evanesce in the medium. Since the nature of the Green’s 
function changes at the boundary z = 0 and z = L, the energy of the evanescent 
wave would return back into the media pertinent to the appropriate region of the 
Green’s function. For example, at z < 0 no energy from the evanescing waves from 
reflection would go into the medium, nor would any energy from the evanescing wave 
from transmission contribute into the air. 

From this point on, we will ignore the effects of the evanescent waves because 
we will assume the transmitter and receiver are sufficiently far away from this region. 
But,we can see that in a half-space z <0, the intensity G*G = 

becomes smaller as z gets increasingly negative. We expect therefore that the evanes¬ 
cent contribution will be small, and we should be able to neglect such effects. 

F. REDUCTION OF THE INTEGRAL 

The motivation behind the coordinate transformations of Equations (IV.47), 
(IV.48),(IV.51),(IV.52) is twofold. First, they will allow us to reduce the double 
integral of Equation V.12 into a single integral, which is more tractible. But we 
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also do this transformation in conformity to Brekhovskikh’s model. This approach 
will allow us to set up the contour integral which we will present in this section and 
perhaps be able to solve in the future. Since the method for which this transformation 
is the same for all the Green’s functions, we will demonstrate it for the case of the 
reflected Green’s function in this section. 

The symmetry in the (j) direction implies a switch to polar coordinates: 

X = pcos(j)' y = psm(j)' (V.27) 


Substituting this equation into (V.12), yields 


G{p,z,ro) = 


>0 Jo 


( — ) cos6*sin6»d6»d(^ 

271” 


Since 


27rJo(/cpsin 9 ) 


We are left with a Bessel Transform: 



^ikpcos(4,-4,') 


2 

27rJo(fcpsin0) cos6*sin6*d0 (V.28) 

Because the Hankel functions are a linear combination of the Bessel series of 
the hrst and second order, Brekhovskikh [Ref. 2] writes the Bessel Jo in the following 
way: 


G{p,z,ro) = 


G{9,z,ro) 


k 


Jo{kpsm9) = - Hq \k p sin 9) + HQ^\kp sin 9) 


r(2)/ 


(V.29) 


and, using the Hankel function property 


i/f (e—x) = 


Equation V.28 becomes 


k 


2 PTj—ioo 


1 r 


G{p,z,ro) = ^l G{9,z,ro)^ H^^\kpsin 9) + {kpsin 9) 

ZTT Jq 2 1 


r(2)/ 


cos 9 sin 9d9 
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k‘^ 

47r 



•200 ioo 

G{9, z, ro) cos 6* sin sin 9)d9 + / G(6*, z, Vq) cos 6* sin 0Ho^^(kpsin 9)dd 

Jo 


k^ 

djT 

The reflection Green’s function is therefore 

G{p,z,ro) = — / G{6, z,ro) cos 6 sin 6HQ^\kpsin 6)d6 


G{6, z, ro) cos 6 * sin 6H^\kpsin 6)d6 + 




G(0,z,ro) cos6*sin6*Ho^^(kpsin6*)d6* 


flvr 


m 


^ikzQ cos 0 ! - 

^_ ^—ikzy 1—sin^ 9 

2k cos 9 


cos 9 sin 9Hq^'' {kp sin 9)d9 


_ /■§-*“ (1 - n^) + (n^ - 

flyr (\/n 2 — sin^ 0 + cos 0)2 — (\/n 2 — sin^ 9 — cos 0 ) 2 gi 2 fcL\/n 2 -sin 2 e 

(V.30) 

pikzQ cos ^ ! - 

X —- ^ cos 9 sin 02/^^^ {kp sin 0)d0 

G. ANALYSIS 

Performing the contour integration of Equation V.30 requires very advanced 
mathematical concepts. One has to approximate the integral by assuming certain 
conditions of the argument of the integrand — such as wave number being very large. 
By using the Method of Steepest Decent where one has to And the saddle points and 
determine the right saddle path, an analytical approximate solution can result from 
Equation V.30. We will not perform this integral in this thesis, but we will analyze 
the integral. 

First, the limits of integration tell us that the summation goes from all the 
allowed real angles to the complex angles (which is the region of evanescent waves). 
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Let us examine the latter half of the integral. We saw earlier that Equation V.3 is 
the free space Green’s function. Analogously 


47r 



^ikzQ cos 6 ! - 

—- ^ cos Q sin eH^Q\kp sin 6»)d6» 

cos yJ 


(V.31) 


would give us the free space Green’s function (The only difference being that it is left 
going from the negative sign given as stunning idea is that the 

reflected wave is spherical in nature after performing the Hankel Transform. This is 
the idea behind Huygen wavelets as shown in Figure 11. From the figure you can see 
that the returning waves are spherical, and they superpose with each other to create 
a wave front. The singularities of Equation V.30 are the trapped modes of either 
the TE or TM waves in the media (G. Scandrett, personal communication, December 
2006), [Ref. 32]. This approach can be better seen when converting the Gartesian 
form into the polar form instead of the spherical form (See Appendix G). The branch 
cuts show the forbidden regions, and these are the regions of evanescing waves which 
Brekhovskikh describes in meticulous detail, and whch require a contour integration 
path that will account for these forbidden regions [Ref. 2]. 



® W. Fendt 1998 

© Prof. T. Mzoughi 1998 Angle of Refraction: 14.5 ” 


Figure 11. Huygen Wavelets. Figure was taken from [Ref. 35] 


On the other side of the media {z > L), we can also use Huygen Wavelets. 
From its Gartesian form, we can demonstrate see this effect in a MATLAB simula- 
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tion(see Appendix F). In Figure 12, the observer is on the plane z = 3m, and we can 
see these wavelets. 


Transmission v=500GHz, Observation point at 2*3 



Figure 12. Huygen Wavelets. Observation a,t z = 3m 

It maybe benehcial to perform the integral in Equation V.30 to obtain the 
Green’s function in the spatial domain because it shows us the effects of all these 
spherical waves. However, this thesis will focus on an imaging algorithm, and we will 
work in the spatial frequency domain since we already have an analytical form. We 
will simulate the data and perform the transforms numerically. 

H. MATERIAL ANALYSIS 

Listed in Table II are some common barrier materials. We will examine the 
reflection from the wall. As expected, the type of material will greatly affect the 
amount of reflection from the wall. And there are frequency windows at which we get 
minimal reflection. Figures 13, 14, and 15 show the reflectivity of materials listed in 
Table II, with barrier width of half a meter. 
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Frequency vs Reflection from the dielectric 



Figure 13. Reflection (Cloth, Balsa Wood, and Wooden Door) 

The MATLAB code used to generate these results uses the reflection coefficient 
calculated in this thesis. The code is modeled for a plane wave parallel to the surface, 
and its interaction within the media. 

We see from Figure 13 that the reflection from cloth and wood are very small 
at THz frequencies of current interest and so the material becomes almost transparent 
to these waves. 

As the index of refraction increases, the reflections coming back from the 
medium become more signihcant as shown in Figure 14. But there are key frequencies 
where reflection can be minimized. For the case of Dry Concrete, a common material 
used in construction of most walls, the material essentially becomes transparent at 
500 GHz. But a potential problem relates to the target you trying to image. At 
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Frequency vs Reflection from the dielectric 



Figure 14. Reflection (Plexiglas, Paper, Dry Concrete) 

500 GHz, plane waves will penetrate dry concrete with little loss (at one meter in 
depth). If you want to get information on people, this frequency would be great. If 
you look at Figure 15, muscle is not very reflective at this frequency, but blood is as 
high as 90% reflective. Thus, gathering data scattered from people at this frequency 
would be “ideal”. In general, a balance must be achieved between wave penetration 
of the dielectric and return information from the object—unless the object is metallic 
(in which case the only concern then would be determining the best penetrating 
frequency). 
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Frequency vs Reflection from the dielectric 



Figure 15. Reflection (Diamond, Muscle, Blood) 


Index of Refraction 

Cloth 

1.09 

Balsa Wood 

1.14 

Wooden Door 

1.41 

Plexiglas 

1.51 

Paper 

1.73 - 2.00 

Dry Concrete 

2.12 

Diamond 

2.41 

Muscle at 37 °C 

7.00 

Blood at 37 °C 

7.62 


Table II. Typical index of refraction for common items 
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VI. INTEGRAL EQUATIONS 


A. WAVE EQUATION 

The Electric field obeys the time independent Helmholtz equation. If we allow 
the index of refraction to vary so slowly with position that it is essentially constant 
over distances on the order of a wavelength, then we do not have to concern ourselves 
with the coupling of Maxwell’s Equations [Ref. 25], and thus we may examine its 
scalar form: 

V^E{r) + enl,^{r)E{r) = 0 (VI.l) 

We can re-write this equation as 

V^E(r) + k‘^E{r) = - k^{nobjir)^ - 1) E{r) (VI.2) 

'-V-" 

u 

where U denotes the “scattering potential” of the medium. 

V^E{r) + UE(r) = -U{r)E{r) (VI.3) 

The the Green’s function of the Helmholtz equation will be a solution of 

(V^ + U)G'(r - ro) =-h(r - ro) (VI.4) 

B. LIPPMAN-SCHWINGER EQUATION 

We see from Figure 16 that the Green’s function is divided into three regions 
with solutions discussed in the previous chapter. We also know that the Green’s func¬ 
tion is symmetric so that the response at r due to a concentrated source at Tq is the 
same as the response at Vq due to a concentrated source at r. We will use these ideas 
to derive the Lippman-Schwinger’s Equation and apply the first Born Approximation 
in its solution. In Equation VI.3, E{r) is the total field. If a transmitter is in region 
HI and it emanates a planar wave Eindx) ~ into the media, then this incident 

wave will be planar inside the medium (region H), and comes out as a planar wave on 


59 



the other side into region I as = tAe where t is the transmittance. This 

held satishes the homogeneous Helmholtz equation 

V^EUr) + k^EUr) = 0 (VL5) 

E[r) can be expressed as the sum of incident held and the scattered held so that 

E(r) = + BUr) 

So where Egdr) is the held scattered from the target. Equation VI.3 can be therefore 
simplihed to 

V^E,,{r) + eE,,{r) = -U{r)E{r) (VI.6) 

Multiplying Equation VI.4 by Esc{r) and Equation VI.6 by G, and subtracting 
the two yields 

E^Jr)V^G{r - Vf,) - G{r - j'o)V^£sc(r} ^ U{r)E{r)G{r - ro) - E^Jr)S{r - ro) 

(VI.7) 

Using the fact that the Green’s function is symmetric, we integrate with respect 
to ro throughout the entire Volume Vr and apply Green’s Second Theorem to arrive 
at 


Escir)= U{ro)E{ro)G{ro - r)d^ro 
Jvr 


(VI.8) 


'Sr 


(E(ro)V°G(ro - r) - G{r, - r)V°E,,(ro)) • dS 


R 


The Sommerfeld Radiation condition guarantees us that the helds go to 0 at inhnity, 
and therefore, as the surface area increases with the radius R, the area integral 
vanishes. Furthermore, since the Scattering potential is 0 everywhere outside the 
scattering object, the volume integral of the entire volume Vr is then shrunk down 
to just the volume of the scattering object. Therefore, the integral reduces to 


Escir) = f U{ro)E{ro)G{ro - r)d^ro (VI.9) 

Jvo 
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Figure 16. Illustration of a Scattering Medium occupying Volume Vr bounded by a 
closed Surface Sr 

This is the Lippman-Schwinger Equation. The Green’s function is different for each 
region, depending on where the observation point is located, particularly along the 
z-axis: 




G{r - ro) = < 


Giir - ro) 
Giiir - ro) 
Gni{r - ro) 


and Gj, Gjj and G/// are the respective Inverse 
(IV.47), (IV.3), (IV.51), and (IV.52). 


for z < 0 

for 0 < z < L (VI. 10) 

for L < z 

Fourier Transforms of Equations 
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C. BORN APPROXIMATION 


Substituting Equation VI.9 into Equation Byields 


E = + / U{ro)E{ro)G{ro - r)dVo 

JVo 

If we formally define the operator 

{AE){r)= j G(ro - r)E(ro)dVo 

Jvo 

E = E'-^^ + 

Esc 

= E',„, + A(E',„, + AE) 

= + AEI^ + A\EI^ + AE) 


(VI.ll) 


Then the Born approximation of Equation VI.ll becomes 

E = E[^^+ [ V(ro)ELMG(ro-r)dVo (VI.12) 

JVo 

Esc= [ V(ro)E'„,(ro)G(ro - r)dVo (VI.13) 

JVo 


D. THE MODELING SCHEME 

In this thesis, we will only account for single scatter events, assuming the scat¬ 
tering from other points to be signihcantly smaller and so will not add any meaningful 
contributions. We shall now the Scattering Potential. For point scatterers, we can 
easily model this as a series of Dirac Delta Functions so that 


U{r) = '^Ai6{r - Vi) (VI.14) 

i 

where Ai are (complex-valued) scatterer amplitudes. Substituting Equation VI. 14 
into Equation VI.13 and integrating, we get a summation of Green’s functions. Since 
the receiver will be in region III, we use the Transmitted Green’s function in region 


III. 


Esc{r) = AAte-^’^^'Gniir - r,) (VI.15) 

i 
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VII. 


ALGORITHM 


A. SPATIAL DOMAIN VS FREQUENCY 

The Green’s function forms a kernel which operates on the Scattering Potential 
and the Electric held. It is also shift invariant, giving us the special characteristic of 
the convolution integral. Let us revisit and focus on the volume integral of Equation 


VI.8 


Escir) = f U{ro)E{ro)G{ro - r)d^ro (VII.1) 

JVr 


When R becomes inhnitely large, the volume integral extends over all space. 


CX) 

U (xo, yo, zq)E{xq, yo, ^;o)G'(x-xo, y-yo, z- 2 ;o)dxodyodzo (VII.2) 

— oo 

It is the scattering potential that limits this integration to the scattering volume. The 
hrst Born approximation of this integral yields 


Esc{x,y,z) = 


U (xo, yo, zq)E[^^{xo, yo, zo)G{x - xq, y - yo, - 2 ;o)dxodyodzo 

(VII.3) 

Since we know the kernel in terms of its spatial frequencies, we perform the 
convolution integral in the Fourier domain. We reduce the 3D Fourier Transform to 
a 2D Fourier Transform by holding z constant. The idea is that the receiver will 
be placed at a hxed z measurement. The 2D Fourier Transform of this convolution 
integral will result in the 2D Fourier Transform of the Green’s function, and a 2D 
Fourier Transform of the product of U and which is U convolved with E^^^ 

in the frequency regime. Since our scattering model is a summation of Dirac Delta 
Functions, the 2D Fourier Transform of Equation VI. 15 with respect to x and y is 
simplihed to 


E,c{x,y,z) = 
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(VII.4) 


Tni{x - Xi,y - yi)e 
The scattered field in the spatial frequency domain for a given z is then 


/ /*oo _ 

/ AAte-^^^'GuAx -x„y- 

' — OO 


E{f., fy) = Y. AAte-^’^^Giuif., 

i 

The advantage of using the frequency domain is that we deal with only prod¬ 
ucts in MATLAB, and therefore, no fancy coding will be required (simple element-to- 
element matrix multiplications are enough). We will also let the computer calculate 
the 2D Inverse Fourier Transform to the spatial regime and compare the results. 

B. MATLAB CODE IMPLEMENTATION 

Equation VII.5 is the equation we shall use to create our simulated scattered 
field data. The total field at the transmitter will be composed of the transmitted 
field of the device, the reflected fields from the wall, the scattered field, and the 
lateral (evanescing) field. We assume that the evanescing waves will be minimal, and 
that we are purposely staying away from the evanescing regions. Therefore, the only 
factors we account for are the transmitted, reflected, and the scattered fields. For 
our simulation, we neglect the details of the transmitted fields and reflected fields 
as these will be known measurements in the real world. For image processing, we 
are interested in the scattered field information, and whether Tikhonov and/or SVD 
methods will work in mitigating noise corruption. 

For simulation purpose, we reduce the propagation number so that the effects 
of the point scatterers will be amplified. The limitation of this model is that it is 
computation intensive, and we are limited to the sampling frequency and bandwidth. 
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Point Scatterers 



Figure 17. Amplitude of field measurements due to 5 Scattering Points 

Figure D displays the noise-free field measurements, collected from 5 scattering 
points measured at z = 2m. To this data, we add Gaussian noise, and we will apply 
the method of Tikhonov and SVD hltering discussed in Chapter 2. 

The code is presented in Appendix D. The idea is to do all the work in the 
spatial frequency and use a 2-D inverse FFT to bring the results back to the spatial 
domain. 

C. SIGNAL TO NOISE RATIO 

The additive noise is Gaussian noise and applied in the frequency domain. 
The SNR is calculated by summing up the intensities of each element of the matrix 
for the signal and the noise, and taking the ratio of the two. This provides a way to 
gauge our data. We increased the Gaussian noise by multiplying it by a noise factor 
which gives us the desired SNR value. 
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Noisy Data with SNR of -30.00 



Figure 18. 5 Scattering Points under Gaussian Noise SNR=—30dB 

D. SPATIAL FREQUENCY DOMAIN 

Consider the case of a two-point scatterer. For one scatterer, the frequency 
domain follows Figure 19. The second point scatterer, looks similar to the hrst one 
(see Figure 20). The phase shift is hard to observe. 

However, when we add the two, we can see the effects of phase interference, 
(see Figure 21, and upon Fourier inversion, we get what we are supposed to get - two 
scattering points (see Figure 22). 

E. TIKHONOV REGULARIZATION 

From Figure 18, the noise has almost completely masked the signal. Now we 
apply Tikhonov regularization, and the we will choose the regularization parameter by 
trial and error. We know that as the regularization parameter decreases, the linear 
compact operator discussed in Chapter 2 will head towards the Least Squares 
solution to the problem — and thereby increase the condition number. Thus, we will 
see that the recovered image will become more unstable (see Figure 23). 
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Single Point Scatterer in Frequency spectrum (Transmission) 



Figure 19. Single Point Scatterer (Spatial Frequency Domain) 

This figure shows the ill-conditioned nature of the problem. When the reg¬ 
ularization parameter is increased to a = 10, we see how well the noise can be 
mitigated (see Figure 24). There is a trade-off between a large regularization param¬ 
eter and a small one. A larger regularization parameter will dominate the inverse 
operator giving you optimization problems, while a small one will give you a least 
squares solution but a diverging condition number. 

F. TRUNCATION SVD 

We now apply the SVD-truncation method. At first, the singular values were 
truncated at the very large value a = 10000, and we see the ill-conditioned behavior 
(see Figure 25). 

The idea is to compare the singular values, and try to locate where the rapid 
changes of singular values occur and truncate at that point. One parameter was 
determined to be at a = 1, and the image is pulled out of the noise (see Figure 26). 
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Figure 20. Second Point Scatterer (Spatial Frequency Domain) 

One can see that the SVD method appears to do a visually better job in 
retrieving the signal from the noise, but this may not be generally true. 


Second Point Scatterer in Frequency spectrum (Transmission) 
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Figure 21. Two Scattering Points (Spatial Frequency Domain) 
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Figure 22. Two Scattering Points (Spatial Domain) 


Two Point Scatterer in spatial domain (Transmission) 
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Tikhonov Regularization alpha=0.003 and SNR=-30.00 


Figure 23. Tikhonov Regularization a = 0.003 at SNR=—30 dB 


Tikhonov Reguiarization aipha=10.000 and SNR=-30.00 


Figure 24. Tikhonov Regularization a = 10 at SNR=—30 dB 
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SVD Truncation w/ Singular values truncated at 10000.00 
SNR=-30.00 


SVD Truncation w/ Singular values truncated at 1.00 
SNR=-30.00 


Figure 26. TSVD a = 1 at SNR=- 


Figure 25. TSVD a = 10000 at SNR 
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VIII. 


CONCLUSION 


Our model has assumed linearly polarizable media with unchanging permit¬ 
tivity. For non-linear materials, the polarization vectors will have to be accounted 
for separately. Though tedious, this extension would provide us with a more general 
model. In addition, the material we considered is not magnetizable. These assump¬ 
tions are generally good, but cannot be true for all Through-the-Wall cases. It is 
important to note that when the media is both non-linear and magnetic, the Maxwell 
equations will be coupled and the scalar theory has to be readjusted to ht this model. 
A complementary examination of common wall materials will be required, and some 
of the future work is to see how well the model works in “real-life” applications, and 
how we can adjust the model to ht specihc problems. We have demonstrated in this 
thesis that the scattered held can be modeled as simple summation of Green’s func¬ 
tions because we applied the the weak scatterer model. When the Born approximation 
is not valid, we may have to consider other approximations, and we may also want 
to include multiple scattering in which case the Lippman-Schwinger equation may be 
extended to include for scatterers from nearby points. 

The Rytov Approximation assumes that the wave perturbation is caused by 
objects through a phase change in the incident wave along the propagation path. 
It is derived similarly to the way we derived the Born approximation, except the 
held is comprised of the phases as functions of position. It is said that the Rytov 
approximation is a more accurate approximation than the standard Born model. In 
future work, one may take the Lippman-Schwinger Equation and apply Rytov and 
compare to the Born model. It may be worth applying this to “real” world problems 
along with the Born model for future work. 

We discussed non-conductive dielectrics. For conductive media, the attenua¬ 
tion maybe too great for examining the scattered held. The attenuation of a signal 
through the media and back maybe masked by the Gaussian noise and the data may 
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be too noisy for any reasonable image reconstruction algorithm to work. For direct 
application to real world scenarios, the conductive case may be unlikely to be en¬ 
countered. For imaging behind walls, most construction materials are not completely 
made of conductive materials. However, steel bars may pose a problem with back 
scatter inside the media, which we would have to deal with separately. 

We also ignored dispersive effects that need to be examined for various mate¬ 
rials for which a priori knowledge of the material characteristics must be known in 
order to make the appropriate approximations for the theory. It will be of interest in 
future work to make appropriate approximations and account for dispersion by using 
the Method of Stationary Phase. Such approximate analytical solutions can provide 
some understanding about how the waves behave in transit through the medium, but 
also can be much more efficient and accurate in the image reconstruction stage. 

The Tikhonov regularization parameter was found by trial and error. There 
are much more effective ways to determine the regularization parameter such as the 
L-curve method [Ref. 30]. The Bayesian statistical approach [Ref. 30] to the problem 
can give us better image reconstruction and might be fruitful for future research. 

Our code implementation has some limitations. Since MATLAB can be very 
memory hungry, the bandwidth examined could not be very wide, nor can one increase 
the sampling frequencies very high without sacrihcing time (and patience). A more 
efficient computation scheme will be required. In this thesis, I have used short band- 
widths with sampling frequencies at Nyquist frequencies to avoid aliasing. A dual 
core 2Ghz with 6GB of RAM would take 5 — 15 minutes to do one full computations 
for a reasonable frequency bandwidth of 20m“^ for and fy. This poses a problem in 
computation and analysis on systems requiring higher bandwidths. Efficient coding 
and higher computing power will be required for analysis of larger bandwidths. 
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APPENDIX A. THE RIEMANN LEBESGUE 

LEMMA 


lim / K(A, s)C'sin(ns)(is = 0 


(A.l) 


Through integration by parts, we get the following and therefore conclude that as 
n —> oo, the integral will then goes to 0. 


lim — 

n—>oo 


Ks) cos{ns) 1 
T o 


n 


a n 


dKi\,s) , ^ 

, ---cosfnsjas = 0 

' os 


(A.2) 
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APPENDIX B. CURRENT IMAGING 
METHODS (ANL) 


At Argonne National Laboratory, A receiver and transmitter via waveguides 
are used for imaging. The wave guides are shown in the following diagram. The 
conhguration is set to send electromagnetic waves from 10 Ghz to 200 Ghz, but it 
also has the ability to go to higher frequencies if necessary. Both the transmitter and 
receiver are of the same type and design, and they are adjacent to each other while 
in radar mode (i.e. data received from scattered Helds). In the radar mode, both the 
transmitter and receiver are moved simultaneously in space in sequence. 



Figure 27. Waveguide for transmitter. Separate waveguide for receiver not presented 
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APPENDIX C. SINGULAR VALUE 
DECOMPOSITION 


Let A be a linear operator that maps X ^ Y, and to be more specific, let 
us say that X C and Y C So that A will have dimension m x n. In this 
section we will discuss (n x n) and AA"^ {m x m). 

The eigen equation for A^A is 

A^ Axi = XiXi (C.l) 

where \i is the eigenvalue corresponding to the vector Xi. Because A^A is symmetric, 
the eigenvectors are orthogonal [Ref. 29]. Similarly, the eigen equation for AA^ is 

AA'^r/j. = (C.2) 

where xl^j is the eigenvalue corresponding to Uj. AA^ is also a symmetric matrix and 
therefore the eigenvectors y^are orthogonal. 

Consider the eigenvector Xi of A"'"A, with nonzero eigenvalue A*. Then Ax^ ^ 
0 because Axi is an eigenvector of AA"*": 

{AA^)Axi = A{A^A)xi = A,(A®i) (C.3) 

The same is true for an eigenvector y of AA^ with nonzero eigenvalue since j^y is 
an eigenvector of A'^A. 

It follows that the non-zero eigenvalues of A'^A must also be eigenvalues of 


79 



AA^ and vice versa. If the rank of A"^A is r, then this means that 


Ai = -01 

A 2 = 02 

\r = ll^r 

Ar+1 0 0J-+1 0 


Therefore 


An 0 0m 0 


X = A^y y = Ax 

And normalizing the eigenvectors, we arrive to the following conclusion: 


X = 


_ A^y 


Ax 


jA'j/ii ^ \\Ax\\ 

The singular values Ok of A are dehned by 

\\Axk\\ = = cTfc = a/A^ = \f¥k 

Then equation C.5 becomes 


(C.4) 


(C.5) 


(C.6) 


A^y,^ = UkXk Axk = akyk (C-7) 

Since the eigenvalues of A^ = 0 and 0^ = 0 for fc > r, this means that 
Axk = 0 forfc = r + l,...n 
A^y^ = 0 for fc = r + 1,...m 

If we take Ax^ = (Tky^., and multiply equation C.7 by x^ on both sides, and 
use equation C.8, we see that 

r 

A = '^a^yixJ (C.8) 

i=l 
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and similarly, 


r 

= aiX,yJ (C.9) 

i=l 

The orthonormal vectors x are called right singluar vectors and y are called 
left singular vector. In matrix form, we can write equation C.8 as 


or 

A = YSX^ (C.ll) 

where /S' is a diagonal matrix of the singular values, and Y and X are composed of 
the eigenvectors corresponding to those singular values. 

The inverse is easily calculated by the following: 

A-^ = {YSX^)-^ 

Since X"^ and Y are orthogonal matrices due to the fact that the eigenvectors are 
orthogonal, we see the following: 

A-^ = (X^)-^S-^Y-^ 

A-^ = XS-^Y'^ (C.12) 

Combined with the Least Square Method, this result can be used to get us a solution 
to the system of equation. As you can see the inverse mapping provides us with ^ 
for its diagonal matrix, and this becomes a problem as the singular values get smaller 
and smaller. Since the non-zero eigenvalues of a Hilbert Space operator are isolated 
and form a sequence that converges to 0 [Ref. 28], the norm of inverse mapping 
A~^ whose singular values of will increase indefinitely. Therefore, regularization 
methods need to be introduced to cure this effect. 
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APPENDIX D. MATLAB CODE FOR 
SCATTERING POINTS 


clear; 


"/oParameters 
nu=500e9; 
lambda=3.OeS/nu; 
yok=2*pi/l ambda; 
k=3000; 
n=2.4; 

L=0.5; 
z=2; 

A=l; 

signal=0; 
noise=0; 


% Frequency of the wave propagation. 

°/o In terms of wave length. 

% Magnitude of the propagation vector. 

% Material index of refraction. 

% Thickness of the wall. 

°/o Measurement plane at z= constant. Adjust accordingly. 
°/o depending on position of the receiver. 

% Amplitude of the incident plane wave. 


% Defining Coordinate Systems in the spatial frequency domain. 
ii=20; °/o spatial frequencies for meshgrid input. 

fb=ii; % frequency band. 

y„fn=fb; 


fn=2*fb; y Nyquist rate. 

j=l/fn; y Separation must be at least this distance to prevent aliasing, 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii); y Creating coordinate system. 


y Defining the Green’s function. 

x0=0;y0=0;z0=-3; y Set this position accordingly to the problem. 

sN=sqrt(k~2-(2*pi*fx) . "2-(2*pi*fy) . "2) ; °/„ Look in chapter 3 for reference. 

sNm=sqrt ((n*k) ~2-(2*pi*fx) . ~2-(2*pi*fy) . ■'2) ; y Look in chapter 3 for reference. 
num=4*sN.*sNm*exp(-i*sN*L); y Numerator of the equation. 

y Denominator of the equation. 

den=(sMm+sN).~2.*exp(-i*sNm*L)-(sNm-sN)."2.*exp(i*sNm*L); 

X0=(exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0).*exp(-i*sN*zO).*exp(i*sN*z))./(2*sN); 


G=(num./den).*X0; y Transmission Green’s function. 
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T/X/X/XyXtyXyXyXTLTL°L synthetic Data Generation 


0/0 /o /o/o 


% transmittance for a plane wave through a given media 
t=4*n*exp(-i*k*L)./((n+1)~2*exp(-i*n*k*L)-(n-1)~2*exp(i*n*k*L)); 


% Scattering Points in the spatial frequency domain. 

xl=-l;yl=-l;zl=-l; 

x2=-l;y2=l;z2=-3; 

x3=0;y3=0;z3=-5; 

x4=l; y4=l;z4=-1.3; 

x5=l; y5=-l;z5=-2.2; 

% Incident field after the penetration of the media 
El=t*A*exp(i*k*zl); 

E2=t*A*exp(i*k*z2); 

E3=t*A*exp(i*k*z3); 

E4=t*A*exp(i*k*z4) ; 

E5=t*A*exp(i*k*z5) ; 


% Positions of the scatter points in frequency domain for the given 
% measurement plane at z=constant 

% Since the scatter points are in terms of dirac delta functions, 

% the field measurement would become the summation of the Green’s function 
% and the incident wave after performing the Born Approximation 


Xl=(exp(-i*2*pi*fy*xl).*exp(-i*2*pi*fx*yl).*exp(-i*sN*zl).*exp(i*sN*z))./(2*sN); 
X2=(exp(-i*2*pi*fy*x2).*exp(-i*2*pi*fx*y2).*exp(-i*sN*z2).*exp(i*sN*z))./(2*sN); 
X3=(exp(-i*2*pi*fy*x3).*exp(-i*2*pi*fx*y3).*exp(-i*sN*z3).*exp(i*sN*z))./(2*sN); 
X4=(exp(-i*2*pi*fy*x4).*exp(-i*2*pi*fx*y4).*exp(-i*sN*z4).*exp(i*sN*z))./(2*sN); 
X5=(exp(-i*2*pi*fy*x5).*exp(-i*2*pi*fx*y5).*exp(-i*sN*z5).*exp(i*sN*z))./(2*sN); 

% Generating Data 

h=(num./den).*(E1*X1+E2*X2+E3*X3+E4*X4+E5*X5); 


"/oSetting matrix dimension 
[M,N]=size(h); 

h=h(l:M, 1:N) ; "/oChange matrix dimension if warranted. 

"/oConverting bin numbers to real axis numbers in spatial domain (meters) . 
"/oRange has to be up to Nyquist frequency. 
ax=linspace(-fn/2,fn/2,length(h)); 
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n=(randii(M,N)+i*randn(M,N)) ; “/onoise in the frequency domain. 


% Set the signal-to-noise (in dB) 
snr=-35; 


"/oSignal to Noise ratio 

SIG=abs(h.*conj(h)); 

N0I=abs(n.*conj(n)); 
signal=sum(sum(SIG)); 
noise=sum(sum(NOI)); 
factor=sqrt((signal/noise)*10" 
n=n*factor; 

d=h+n; 


-snr/10)); 

“/odata in the frequency domain. 


yoyoyoyoy.yormyoyoyoyoyoyoyoyoyoyoy.y.yoyo End of (noisy) Data Generationyoy;/;/;m%%%y„yoyoyoy„yoyoy.yo 


figure(1) % data corrupted with noise in spatial domain. 

dl=ifft2(d); % data in the spatial domain. 

dl=ifftshift(dl); 

imagesc(ax,ax,abs(dl)); 

colormap (1-gray); 

axis([-2 2-22]) 

title(sprintf (’Noisy Data with SNR of yo2.2f ’ , snr)) 

% Tikhonov Regularization 

alpha=3; % Regularization parameter for Tikhonov. 

R=inv(alpha*eye(size(G’*G))+G’*G)*G’; 

figure(2) 

p=R.*d; 

p=ifft2(p); % data in the spatial domain. 

p=ifftshift(p); 

imagesc(ax,ax,abs(p)); 

axis([-2 2-22]) 

colormap (1-gray); 
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title(sprintf (’Tikhonov Regularization alpha=°/o2.3f and SNR=°/o2.2falpha, snr) ) 


% SVD Truncation Method 

alpha=.l; % Parameter for Truncation of singular values, 

[u,s,v]=svd(G); 

% inverse of the singluar value matrix. 
sinv=zeros(N,M); 
for l=l:rank(s) 

if 1/s(1,1)<alpha 

sinv(l,l)=l/s(l,l); 

else 

sinv(l,l)=0; 

end 

end 


e=(v*sinv*u’).*d; 

e=ifft2(e); % data in the spatial domain. 

e=ifftshift(e); 


figure(3) 

imagesc(ax,ax,abs(e)); 
colormap (1-gray); 
axis([-2 2-22]) 

title (sprintf (’TSVD w/ singular values cutoff at yo3.2f \n SNR=°/o2.2falpha, snr)) 


86 



APPENDIX E. MATLAB CODE FOR 
MATERIAL REFLECTION 


clear; 

"/oCalculating the Reflection for plane waves in various surface 
"/oEquations used from Chapter 3 

"/oCreating various index of refraction for different materials 

"/oCloth, Balsa Wood, Wooden Door, Plexiglas, Paper, Dry Concrete, Diamond, 

"/oMuscle, Blood 

nm=[1.09 1.14 1.41 1.51 1.73 2.12 2.41 7.00 7.62]; 

"/oThickness of the material 
L=0.5; 

c=3.0e8; °/oSpeed of light 


"/oCloth, Balsa Wood, Wooden Door 
figure(1) 
for j=l:3 

nu=linspace(20e9,lel2,200) ; 
k=2*pi*nu./c; 

num=((l+nm(j))*(l-nm(j))+(nm(j)+l)*(nm(j)-1)*exp(i*2*k*nm(j)*L)); 
den=((nm(j)+l)*(nm(j)+l)+(nm(j)-1)*(l-nm(j))*exp(i*2*k*nm(j)*L)); 
r=num./den; 

R=r.*conj(r); 
plot(nu,R); 

text (nu(30*( j-l)+8) ,R(30* (j-l)+8) , sprintf (’n= °/o2.2f ’ ,nm(j))) ; 
hold on 


end 

hold off 

xlabel(’frequency (\nu) Hz’); 
ylabel(’Reflection’); 

title(’Frequency vs Reflection from the dielectric’) 
text (lell,. 115, sprintf (’L = °/ol.2f m’,L)); 
“/oaxisCEO.lell lOell 0 1]); 
clear figure 


"/oPlexiglass,Paper, Dry Concrete 
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figure(2) 
for j=l:3 

nu=linspace(20e9,lel2,100) ; 
k=2*pi*nu./c; 

nuiii=( (l+ninC j+3) ) * (l-nm( j+3) ) + (nm( j+3)+l) * (nin( j+3) -1) *exp(i*2*k*iim( j+3) *L) ) ; 
den=( (nm( j+3)+l) * (niii( j+3)+l) + (nm( j+3) -1) * (l-nm( j+3) ) +exp(i+2+k*nm( j+3) *L) ) ; 
r=num./den; 

R=r.+conj(r); 
plot(nu,R); 

text (nu(30*(j-l)+10) ,R(30* (j-l)+10), sprintf (’n= °/„2.2f ’ ,nm(j+3))) ; 
hold on 


end 

hold off 

xlabel(’frequency (\nu) Hz’); 
ylabel(’Reflection’); 

text(lell4,sprintf(’L = °/ol.2f m’,L)); 

“/oaxisCEO.lell lOell 0 1]); 

title(’Frequency vs Reflection from the dielectric’) 


"/oThickness of the material 
L=0.01 


"/oDiamond, Muscle, Blood 
figure(3) 
for j=l:3 

nu=linspace(20e9,lel2,100) ; 
k=2+pi+nu./c; 

num=((l+nm(j+6))+(l-nm(j+6))+(nm(j+6)+l)*(nm(j+6)-1)+exp(i+2+k*nm(j+6)*L)); 
den=((nm(j+6)+l)*(nm(j+6)+l)+(nm(j+6)-1)*(l-nm(j+6))+exp(i+2+k*nm(j+6)*L)); 
r=num./den; 

R=r.+conj(r); 
plot(nu,R); 

text (nu(30*(j-l)+10) ,R(30* (j-l)+10), sprintf (’n= y„2.2f ’ ,nm(j+6))) ; 

°/otext (nu(6*( j+2) ) ,R( (j+5)) , sprintf ( ’n= %! . 2f ’ ,nm( j+6) ) ) ; 
hold on 


end 



hold off 

xlabel(’frequency (\nu) Hz’); 
ylabel(’Reflection’); 

title(’Frequency vs Reflection from the dielectric’) 
text (lell, . 95, sprintf ( ’ L = °/ol.2f m’,L)); 
axis([0.1ell lOell 01]); 
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APPENDIX F. MATLAB CODE FOR 
TRANSMITTED WAVE FRONT 


clear; 

"/oThis code graphs the transmission of the waves 
"/oObservation of the fields at (0,0,3) 

"/oParameters 

nu=500e9; % Frequency of the wave propagation. 

lambda=3.0e8/nu; % In terms of wave length. 

k=2*pi/lambda; % Magnitude of the propagation vector. 

yok=3000; 

n=2.4; % Material index of refraction. 

L=0.5; % Thickness of the wall. 

z=3; °/o Measurement plane at z= constant. Adjust accordingly. 

% depending on position of the receiver. 

signal=0; 
noise=0; 

% Defining Coordinate Systems in the spatial frequency domain. 
ii=1000; °/o spatial frequencies for meshgrid input. 

j=.89; 

yoj=l/fn; % Separation must be at least this distance to prevent aliasing, 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii); % Creating coordinate system. 

[x,y]=meshgrid(-ii:j:ii,-ii:j:ii); 

% Defining the Green’s function. 

x0=0;y0=0;z0=-3; % Set this position accordingly to the problem. 

sN=sqrt(k~2-(2*pi*fx)."2-(2*pi*fy)."2); % Look in chapter 3 for reference. 

sNm=sqrt((n*k)~2-(2*pi*fx).~2-(2*pi*fy)."2); % Look in chapter 3 for reference. 
num=4*sN.*sNm*exp(-i*sM*L); % Numerator of the equation. 

% Denominator of the equation. 

den=(sNm+sN).~2.*exp(-i*sNm*L)-(sNm-sN)."2.*exp(i*sNm*L); 

X0=(exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0).*exp(-i*sN*zO).*exp(i*sN*z))./(2*sN); 

G=(num./den).*X0; % Transmission Green’s function. 

figure(1) 
imagesc(abs(G)) 
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colormap(1-gray) 

title(’ Transmission \nu=500GHz, Observation point at z=3’) 
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APPENDIX G. GREEN’S FUNGTION POLAR 
GOORDINATE TRANSFORMATION 


If we make the following coordinate transformation: 


f! + fy=p" 


f^ = pcos(j) fy = psm> 


2,2 2 

X + y = r 


X = r cos 6 y = r sin 9 

Then the reflection Green’s function of the Equation IV.47 can be written as 


G{r,9, z,ro) = 27 t / pdp ~ 


X 


Jo(27rpr) k2(n2 - i)(ei2Ly75k^^i7v _ 

2 i|/k“47r2p2 

1 


— 47rp2 -|- _ 4^2p2^2 _ (y/7^2/j.2 _ 47i-2p2 _ y//j2 _q.jj-2p2^2g*2Lfc2 477^p2 

Let 


2ttp nku nkdu 

u = — p= dp = 


nk 


27r 


27r 


G = H - 1) I 


MduJo(nkru) ~ ^ zo)k^/^ 




2/1/2 




X 


fc^(n\/l — + \/l — 'n?v?)‘^ ~ kP‘{n\/l — — \/l — 7^,2^2^2g^2L?^fc^/l 


G = 


X 


kn^i^n? — 1 ) 4 iduJo(nkru) e zo)kVi 


dvr 


a/1 — 


2/J/2 




(n\/l — V? + \/l — 'n?v?)‘^ ~ (?r\/l — v? — \/l — 7 ^ 2 ^ 2 pg^ 2 L^^fc^/T 
let y4 = n\/l — V? and i? = \/l — n ^«2 0 = nkL^/l — v? 
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Where are the zeroes of the denominator of the integrand? 
where does 


{A + Bf 

+ 2AB 
2AB 
A^ + B‘^ 


Denominator vanishes when 


{A - Bfe 
{A^ + B^ - 

e^2e _ 1 ^ 
gi2e _|_ I 

i2 sin 0 
2 cos 0 


i20 

- 2AB)B^^ 
e*e + e*® 

ttan 0 


ttan(nfcL\/l — = 


_ 2n\/l—u‘^y/l—n‘^z 




Look at RHS Numerator is real and positive for 0 < « < 1/n 
it then becomes imaginary and positive for 1/n < -u < 1 
it then becomes real and negative for u > 1 

Denominator goes to zero at u = u* = yj 


Is 


So RHS, 


«* < 1 for n > 1 


+ 1 
2n2 
+ 1 


2n2 


< 1 


< 1 


> 1 


u* > 1/n? 


n^ + 1 ? 1 
2 n m 


n" + 1 > 2 


n 


> 1 


1 * 1 
- < n < 1 
n 

N 

D 


0 < n < 1/n 
1 /n < n < n* 
u* < u < 1 
u > \ 


real,positive 
positive 

imaginary,positive 
positive 
imag,pos 
neg 

real,neg 

neg 


real > 0 
imag > 0 
imag < 0 
real > 0 
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Look at the LHS i tan(nfcL\/l — v?) 


0 < w < 1 LHS imaginary, either positive or negative 


tan(ia;) = ■itanh(a;) 

u> 1 i ■ itanh{nkL\/u‘^ — 1) = — tanh (nfcL\/ — 1) always negative, real 

The conclusion only possible to have a root for 1/n < « < 1 
What about the end points? 

Always a root at u = 1 0 = 0 

For u = 1/n 


itan(nkL\ 1 --) = itanikLVn'^ — 1) 

V 


2ny^l — l/n^i^/n^u^ — 1 — 1 

n? - I - 1 

as u —> 1/n from above 

tan(fcL\/n^ — 1) = 0 

possible to have a root if kL^/n? — 1 = myr Always a root at n = 1. 

Root at n = 1/n if kLyJ'n? — 1 = myr where m is an integer. 

Range of possible roots 1/n < n < 1 
How many roots are there? 

The “real” equation we need to examine is 

^ (IT A-2\ 2n\/l - v?\/'n?v? - 1 

tan(nfcLVl - u^) = -y—-- 

n^ + 1 — 2n^u^ 

To analyze this equation, lets change variables let 

y = nkL\/l — v? 


u^ = l- 


nkL 


Limits: as n ^ 1, y ^ 0. u ^ 1/n, y ^ nkL-^J 1 — = kL\/n'^ — 1 

Call ymax = kL\/n^ 1 , 0 ^ y ^ ymax 


„2 + l_2n2(l-(=)2) 


Ai r,2 _ 1 _ y'^ 

kL\J ^ [kLY 


— (n^ — 1) + 


2l/2 

[kLY 


2y ^J{kLf{n? - 1) - y^ 


[kiy _(n2 - 1) + 


2y^ 

(kLY 
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RHS 


2y\/yLx-y^ 

ymax^+2y^ 


Thus, we see roots of 


tany 


-J/max + 2j/2 


(Note compilation with the aid of Professor James Luscombe) 


This transcendental equation is analogous to the same transcendental equation 
which solves the TE modes for the Dielectric Slabs [Ref. 32]. The intersection of the 
two curves between the left and right hand side of the equation are the values for 
which propagation in the slab is possible. The branch cuts of the integral represents 
the forbidden regions where propagation is not possible and where the waves begin 
to evanesce. Brekhovskikh and Kong work out in detail of these effects [Ref. 2] [Ref. 
32]. The problem with this integral is the nature of the Bessel function. As we 
integrate to infinity, the Bessel function also head towards infinity, and therefore, 
there is no way to close the contour integral. More clever methods are required, and 
Brekhovskikh model provides a way to perform the contour and to implement the 
contour. Nevertheless, we can see that the poles are the trapped modes in the 
media. 
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